{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Think Bayes: Chapter 3\n",
    "\n",
    "This notebook presents example code and exercise solutions for Think Bayes.\n",
    "\n",
    "Copyright 2016 Allen B. Downey\n",
    "\n",
    "MIT License: https://opensource.org/licenses/MIT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "% matplotlib inline\n",
    "\n",
    "import thinkplot\n",
    "from thinkbayes2 import Hist, Pmf, Suite, Cdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Dice problem\n",
    "\n",
    "Suppose I have a box of dice that contains a 4-sided die, a 6-sided\n",
    "die, an 8-sided die, a 12-sided die, and a 20-sided die.\n",
    "\n",
    "Suppose I select a die from the box at random, roll it, and get a 6.\n",
    "What is the probability that I rolled each die?\n",
    "\n",
    "The `Dice` class inherits `Update` and provides `Likelihood`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class Dice(Suite):\n",
    "    def Likelihood(self, data, hypo):\n",
    "        if hypo < data:\n",
    "            return 0\n",
    "        else:\n",
    "            return 1/hypo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the update looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4 0.0\n",
      "6 0.392156862745\n",
      "8 0.294117647059\n",
      "12 0.196078431373\n",
      "20 0.117647058824\n"
     ]
    }
   ],
   "source": [
    "suite = Dice([4, 6, 8, 12, 20])\n",
    "suite.Update(6)\n",
    "suite.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's what it looks like after more data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4 0.0\n",
      "6 0.0\n",
      "8 0.943248453672\n",
      "12 0.0552061280613\n",
      "20 0.0015454182665\n"
     ]
    }
   ],
   "source": [
    "for roll in [6, 8, 7, 7, 5, 4]:\n",
    "        suite.Update(roll)\n",
    "        \n",
    "suite.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The train problem\n",
    "\n",
    "The Train problem has the same likelihood as the Dice problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Train(Suite):\n",
    "    def Likelihood(self, data, hypo):\n",
    "        if hypo < data:\n",
    "            return 0\n",
    "        else:\n",
    "            return 1/hypo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But there are many more hypotheses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0028222671142652746"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hypos = xrange(1, 1001)\n",
    "suite = Train(hypos)\n",
    "suite.Update(60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the posterior looks like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEACAYAAACtVTGuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4HNWd7vHvT2rtkiXLljfJm7yCAxhMjBmSIEIIBgIm\nzAbJzULu3CELuXOTuRPILMGeJZPJfbIRJhfIMJkhkwyZZGDiBHC4BJRAEoxjYwxe5V2yLVu29n07\n949utdqN1F2yutUt1ft5nn5cVV2n+1Qh9OqcU6fKnHOIiIiMJiPVFRARkfSmoBARkZgUFCIiEpOC\nQkREYlJQiIhITAoKERGJyVNQmNl6M9tnZgfM7L5R9nnQzGrMbKeZrfZS1sw+bWZ7zewNM/vS+A5F\nRESSIRBvBzPLAB4CrgdOAtvM7MfOuX0R+9wELHHOLTOzq4CHgXWxyppZFXArcIlzrt/MZib64ERE\nZPy8tCjWAjXOuWPOuT7gCWBD1D4bgMcBnHNbgWIzmx2n7CeALznn+kPlzo77aEREJOG8BEU5UBux\nXhfa5mWfWGWXA+8ys1fM7EUzu3IsFRcRkYkRt+vpApnH757unFtnZm8H/gOoTFJ9RETkAnkJihPA\ngoj1itC26H3mj7BPdoyydcCTAM65bWY2aGYznHPnIj/YzHQzKhGRC+Cc8/JHe1xeup62AUvNbKGZ\nZQN3Apuj9tkMfBjAzNYBzc6503HK/hfw7lCZ5UBWdEgMcc4l5HX0xDkOHT/D4OBgwj5zIl8PPPBA\nyuuQLi+dC50LnYvYr0SK26Jwzg2Y2b3AcwSD5THn3F4zuyf4tnvUOfeMmd1sZgeBDuDuWGVDH/3P\nwD+b2RtAD6GgSZYDR0/z+a89BcCf//FNrFm1MJlfJyIyZXgao3DObQFWRG17JGr9Xq9lQ9v7gA95\nruk4feepX4eXv/jos/znNz4+UV8tIjKp+WZmdl//QKqrMG5VVVWprkLa0LkYpnMxTOciOSzRfVmJ\nZmYuEXX84iPPsn3PsfC6WhQiMpWZGW4CB7OnhJJpeamugojIpOSboJheXHDeerq3pERE0oVvgiI7\nK/O89baO7hTVRERkcvFNUERrau1MdRVERCYF3wRFdE9TY4uCQkTEC98ERbRmtShERDzxbVA0tnak\nugoiIpOCb4OiSV1PIiKe+Dgo1KIQEfHCN0ERPW+iqa0rRTUREZlcfBMU0Rqb1aIQEfHCt0HR1Nap\n2dkiIh74JiiiQ6G/f4COrt4U1UZEZPLwTVCMpFED2iIicfk6KM5pnEJEJC6fB0V7qqsgIpL2fB0U\nDU0KChGReHwTFCNd4HRWQSEiEpdvgmIkZ5vaUl0FEZG05/OgUItCRCQe3wSF4619T2ebOzTpTkQk\nDt8ExUj6+wdobdcjUUVEYvFNUIzWcGho1DiFiEgsvgmK0ZzVXAoRkZgUFBrQFhGJyZdBEQhkhpd1\nGw8Rkdg8BYWZrTezfWZ2wMzuG2WfB82sxsx2mtnqeGXN7AEzqzOzHaHX+vEfjjdl0wvDy5qdLSIS\nW9ygMLMM4CHgRmAVcJeZrYza5yZgiXNuGXAP8LDHsl91zl0Rem1JxAGNKmI0e2ZEUGjSnYhIbF5a\nFGuBGufcMedcH/AEsCFqnw3A4wDOua1AsZnN9lDWxnsAF6JselF4uaFRLQoRkVi8BEU5UBuxXhfa\n5mWfeGXvDXVV/ZOZFXuu9TiVlhSQYcGMam7rpKe3b6K+WkRk0knWYLaXlsK3gErn3GqgHvhqkuoC\nnD+PIjPDKCsdblWcUatCRGRUAQ/7nAAWRKxXhLZF7zN/hH2yRyvrnGuI2P5t4CejVWDjxo3h5aqq\nKqqqqjxUe3RmxuwZ0zh9rhWA0+damT9n+rg+U0Qklaqrq6murk7KZ3sJim3AUjNbCJwC7gTuitpn\nM/Ap4Admtg5ods6dNrOzo5U1sznOufpQ+TuAN0erQGRQJMrsmUVwILhc39CS8M8XEZlI0X9Eb9q0\nKWGfHTconHMDZnYv8BzBrqrHnHN7zeye4NvuUefcM2Z2s5kdBDqAu2OVDX30l0OX0Q4CRwleLTUh\nhloUQ4ZaFiIi8lZeWhSELl1dEbXtkaj1e72WDW3/sPdqJt7smRFBcVaXyIqIjMY3M7Ojbyc+Ry0K\nERFPfBMU0c5rUZxr1XMpRERG4dugKMjLoSAvB4C+/gGaWjtTXCMRkfTkm6CIbDGE5tpFjVOo+0lE\nZCS+CYqR6MonEZH4fB0Uc2YMz84+pRaFiMiIfBkUFrrDSGTX0ylNuhMRGZEvg2LIvFkl4eWTZ5pT\nWBMRkfTlm6AY6eLX8tmRQdGiS2RFREbgm6AYSXFhHvm52QD09PbR2KLHooqIRPNNUIzUWjCzqO4n\njVOIiETzTVBEsoinZcybNfy8JI1TiIi8lS+DIlL57OHnUKhFISLyVr4JitHGqSNbFCfONE1QbURE\nJg/fBEUki+h7KtcYhYhITL4Mikhzy4rDD/g+c66Vvr6BlNZHRCTd+DIoIlsU2VkBykqDt/JwwEnN\n0BYROY8vgyJa5MS7utMapxARieSboIg163rB3NLw8vFTjRNRHRGRScM3QRHLwnkzwsu1CgoRkfMo\nKFCLQkQkFt8ERWTPU+TMbAiOUQxtqm9oobevf8LqJSKS7nwTFLFkZwWYWxaceOeA2lMa0BYRGaKg\nCFH3k4jIyHwZFIa9Zdv8eQoKEZGR+DIoRhLZoqitV1CIiAzxTVC4EZ9xNywyKI6eOJfs6oiITBq+\nCYp45pUVkxXIBKCptZPmts4U10hEJD34JijiPQ87IyODxRUzw+uHa88mu0oiIpOCp6Aws/Vmts/M\nDpjZfaPs86CZ1ZjZTjNb7bWsmf2pmQ2aWWn0e8kSPY9iSGVkUNQpKEREwENQmFkG8BBwI7AKuMvM\nVkbtcxOwxDm3DLgHeNhLWTOrAG4AjiXkaMapcv5wUBypbUhhTURE0oeXFsVaoMY5d8w51wc8AWyI\n2mcD8DiAc24rUGxmsz2U/RrwZ+M8hoSprCgLLx9S15OICOAtKMqB2oj1utA2L/uMWtbMbgNqnXNv\njLHO42aj9D3NnzOdzMzgKWloaqOto3siqyUikpYCSfrcUUYBQm+a5QF/TrDbKW6ZjRs3hperqqqo\nqqoac4XijGUDEAhksmjeDA6Fup0O153lshUVY/4uEZGJVl1dTXV1dVI+20tQnAAWRKxXhLZF7zN/\nhH2yRym7BFgEvG7BP+8rgO1mttY5dya6ApFBkWyV82eGg+LQ8QYFhYhMCtF/RG/atClhn+2l62kb\nsNTMFppZNnAnsDlqn83AhwHMbB3Q7Jw7PVpZ59ybzrk5zrlK59xigl1Sl48UEhMtcpxCVz6JiHho\nUTjnBszsXuA5gsHymHNur5ndE3zbPeqce8bMbjazg0AHcHessiN9DXG6q8Yr3jyKIZGXyB6p05VP\nIiKexiicc1uAFVHbHolav9dr2RH2qfRSj0QZbR4FBJ92l5mZwcDAIPVnW2lp66K4KG/iKicikmZ8\nMzPbq6yszPNaFQeOnU5hbUREUs83QXH+E+5i93KtWDQnvHzgiIJCRPzNN0ExFisqZ4eX9x2pT2FN\nRERST0ExghWLhoPi4PEGBgYGU1gbEZHU8mVQjPSEu0gzSgqZUVIAQG9fP8dO6vkUIuJfvgwKL1Ys\nHh6nUPeTiPiZb4Ii3hPuokV2P+0/qgFtEfEv3wTFWK1YHBEUhxUUIuJfvgkKrzOzhywun0l2VnA+\nYkNTG2ca25JRLRGRtOeboIgUZxoFELyT7MqIcYrdNSeTWCMRkfTly6Dw6m3L54WX36iJvmGuiIg/\n+DIo4s3MHnLJsuHnM71Zc2LM3VciIlOBL4PCq8qKmeRkZwFwrrmD+rOtKa6RiMjE801QXEhjIBDI\n5OIlEeMUBzVOISL+45uguFBvi+h+0jiFiPiRgiKOty0dHtB+88BJjVOIiO/4Jigu9Bd85fyZFObn\nANDc1qn7PomI7/gmKCJ5vOgJgIyMDC5bOT+8vmNPbRJqJCKSvnwZFGN1xUXDQfHa3uMprImIyMTz\nTVCM5Ql30VZHBMW+w/V0dPUkqloiImnPN0ExHiVF+VTOLwNg0Dl27dfVTyLiHwoKjyK7n3bsUfeT\niPiHL4NijD1PAFx+0YLw8mt7j+syWRHxDV8GxYVYvmhW+DLZptZOao6dSXGNREQmhm+CYqxPuIuW\nkZHB2y9ZFF7fuuvIOGskIjI5+CYoEmHdZZXh5VdeP6zuJxHxBd8ERSJ+p1+6vDx8N9n6s60cP9U0\n/g8VEUlzvgmKSMYFjGYD2VkBrrh4eFB7667DiaqSiEja8hQUZrbezPaZ2QEzu2+UfR40sxoz22lm\nq+OVNbO/NrPXzew1M9tiZnNG+tx0s+7SxeHlV17XOIWITH1xg8LMMoCHgBuBVcBdZrYyap+bgCXO\nuWXAPcDDHsp+2Tl3mXPucuBp4IHEHFJ8Y52ZHemKixcQCGQCcOzkOWrr1f0kIlOblxbFWqDGOXfM\nOdcHPAFsiNpnA/A4gHNuK1BsZrNjlXXOtUeULwAGx3UkEyQ/L5srI7qfXvptTQprIyKSfF6CohyI\nvGVqXWibl31iljWzvzWz48AHgC94r/bYJfIKpXesWRZefml7ja5+EpEpLZCkz/XUt+Oc+0vgL0Nj\nF58GNo6038aNw5urqqqoqqoadwXHY82qBeTnZtPZ3cuZxjb2HznNyspJMcQiIlNUdXU11dXVSfls\nL0FxAlgQsV4R2ha9z/wR9sn2UBbg+8AzeAiKdJCdFeDq1ZX8/JV9APzytzUKChFJqeg/ojdt2pSw\nz/bS9bQNWGpmC80sG7gT2By1z2bgwwBmtg5ods6djlXWzJZGlL8d2DuuI4kjsntoHGPZYe+6crj7\n6dc7D9HfPzD+DxURSUNxWxTOuQEzuxd4jmCwPOac22tm9wTfdo86554xs5vN7CDQAdwdq2zoo79k\nZssJDmIfAz6e8KNLolVL5zGjpIBzzR20dXSzY28tayNu8SEiMlV4GqNwzm0BVkRteyRq/V6vZUPb\nf897NdOPmfGuNct46uc7AXj+13sVFCIyJflzZnYi+p6A66++KLy8Y88xGhrbEvK5IiLpxJdBkShz\ny4q5dHkFAA54PjS4LSIylfgmKJI11eGGa4ZbFS+8so+BgUkxb1BExDPfBEWyrH3bIoqL8gBobOlg\nux6TKiJTjIJinAKBTN69dnis/tlfvpnC2oiIJJ5vgiLyCXcJGssOe+87VoWnou86UMexk42J/QIR\nkRTyTVAk06zSIq6KePrdT6pfT2FtREQSyzdBkez79t123aXh5V/+toam1s7kfqGIyATxTVBEutAn\n3MWyYvEcli2cBcDAwCBbXt6d8O8QEUkFXwZFstx63WXh5S0vvUl3T18KayMikhi+DIpEzcyOtu7S\nxcwqLQKgvbOHZ1/SFVAiMvn5MiiSJTMzgztuuDy8vvnFXWpViMik55ugmKin0F23dgUzpxcC0Nre\nxXO/3jMh3ysikiy+CYqJEghkcsd7hlsV//XznfT0qlUhIpOXgiIJ3n3VSmaUFADQ0tbF07/QWIWI\nTF6+CYrzep6SM5YdlpWVye/ecEV4/cnnX6O1vSu5XyoikiS+CYqJdv26lcwrKwagq7uXHz23I8U1\nEhG5MAqKJAkEMvlvt60Lr295eTenGlpSWCMRkQvjy6BI1jyKaGsvWcTKyjlAcLb293766oR8r4hI\nIvkyKCaKmfGRDVeH13+z8xBvHDiRwhqJiIydb4JiouZRRFu+aDbXXLE0vP7of/yS/v6BlNRFRORC\n+CYoUumjt19Nbk4WACcbWvjxi7oNuYhMHgqKCVBaXMAHblkbXv/hlu2caWxLYY1ERLzzT1C45D3h\nzov171jFovKZAPT1D/DwE79IWXeYiMhY+CcoUiwzM4N7/uCd4bl+r++v4/nf7E1pnUREvFBQTKDl\ni2bzvqrhJ+F956nfqAtKRNKeL4MiGU+48+oD71tL+awSAHp6+/jH77+oLigRSWu+DIpUys4KcO8H\nrwtH1Zs1J/nxC7oKSkTSl6egMLP1ZrbPzA6Y2X2j7POgmdWY2U4zWx2vrJl92cz2hvb/TzObNv7D\nGV06/c2+fNFs3h9xK/Lv/fRVDhw9ncIaiYiMLm5QmFkG8BBwI7AKuMvMVkbtcxOwxDm3DLgHeNhD\n2eeAVc651UAN8PmEHNEk8Yc3XcmyhbMAGBwc5Kv/8jztnT0prpWIyFt5aVGsBWqcc8ecc33AE8CG\nqH02AI8DOOe2AsVmNjtWWefc8865wVD5V4CKcR9NDOk2DhAIZPLZj95Afm42AA1NbRqvEJG05CUo\nyoHaiPW60DYv+3gpC/Ax4FkPdUmIVMyjGMms0iI+eVdVeP3VN47yw59tT12FRERGEEjS53r+VWxm\nfwH0Oee+P9o+GzduDC9XVVVRVVU1nrqllatXV3LLtZfw9C/eAOAHz/6WBXNLWXdZZYprJiKTSXV1\nNdXV1Un5bC9BcQJYELFeEdoWvc/8EfbJjlXWzD4K3Ay8O1YFIoPiQkX26EzUbca9+siGqzl+qjF8\nZ9kH/+1F5pYVs3DejBTXTEQmi+g/ojdt2pSwz/bS9bQNWGpmC80sG7gT2By1z2bgwwBmtg5ods6d\njlXWzNYDfwbc5pzz9ShuZmYGf/rRG5gzM3jhV09vH1989FkaWzpSXDMREQ9B4ZwbAO4leJXSbuAJ\n59xeM7vHzP44tM8zwBEzOwg8AnwyVtnQR38TKAT+n5ntMLNvJfbQJpeiglzu+6P15GQH7zJ7tqmd\nv/m/T9PR5esMFZE0YOl+lY2ZuUTU8Uvf3sK2N48CcN8frWftJYvG/ZnJsGPPcf7+21sYHAxeEHbx\nkrn81SduITsrWcNJIjIVmRnOuYT0s2tmdpq54uIFfOqua8Prew6d4uuP/1wPOxKRlPFNUKR7yylS\n1doVfPB9V4XXt+46wte/+4LCQkRSwjdBMdm8/z2ruTXiTrO/2XmIr3/3BQYGBmOUEhFJPAVFmjIz\nPnL71dz8rreFtwXDQt1QIjKxfBMU58+jSF09xsLM+Ngd15wXFr9+7RBffuw5enr7UlgzEfET3wTF\nZDUUFje9czgstu85xqZvPU1bR3cKayYifuHLoEi3mdnxmBn//Xev4XdvuCK8bf+Rev7qwR9ztqk9\nhTUTET/wZVBMRmbGB963lo/dcU14W219E/d95Uk9y0JEkkpBMcnccu0l/K8PXU9GRvA/XXNbJ3/1\nzc388rcHUlwzEZmqfBMULq2ecTc+77xyGV/4xC0U5ucA0N8/wDe++wL/tvmV8IxuEZFE8U1QTDWX\nLC/nS5+9g4rZ08Pbnvr5TjZ966e6maCIJJRvgmIyzcz2am5ZMX//mfez5uKF4W1v1pzkT7/8I17f\nX5fCmonIVOKboIg0ua55ii0/L5v7/8eN/P76NeHjam3v4m++9VO+95OtmpwnIuPmy6CYajIyMrjz\nprfzwKdupbgoDwAHPPn8a3zuK09y9MTZ1FZQRCY1XwbFZJtH4dUly8v5yud+n0uWDz+W/NjJc3zu\nK0/yo+d26D5RInJBfBkUU9n0afk88Mn38dHbf4esQCYAAwOD/PvTr3L/157i4LEzKa6hiEw2vgmK\nKTiWPSoz49brLuUr9/0+yxbOCm8/XNvA/V99km//8CU9OU9EPPNNUPhR+awS/u5PbueD77uKQKh1\n4YAtL+/m03/3BC/9tmZKXg0mIomloJjiMjMzuOOGy/n6/X/A6pXzw9tb2rr4+nd/zue/9hT7Dten\nsIYiku58ExSRfzlP1cHsWOaWFfOXH7+Zz370BqZPyw9vrzl2hr/4xn/xfx77GacaWlJYQxFJV4FU\nV0AmjplxzeVLuHzlfH703HZ++os3wldCvbLrCNt2H+M96y7ijhsuZ+b0whTXVkTShW9aFDIsPy+b\nD2+4mm/+xZ1cc8XS8PaBgUF+9qvdfPJvvs+3f/iSbmEuIoCPgmIyPuEu2WbPmMZnP/Ie/v4z72dl\n5Zzw9oGBQba8PBwYp8+1prCWIpJq6noSli+azd/+zw28vr+OJ57ZRk1orsVQYPzs5d1cffkSNlx3\nGUsjLrcVEX/wZVD4cTA7HjNj9cr5XLaigp376vjBs8OB4Qg+q/vXrx3i4iVzue3dl7Hm4gXhZ2KI\nyNTmy6CQ0ZkZl180n9UrK3h9fx0//vnr7DowfCfaPYdOsefQKcqmF/Heay7m+nUrw/eXEpGpSUEh\nIxpqYaxeOZ8jdWfZ/OLrvLzjUPjBSA1NbXzvp1t54tltrLuskvXvWMVFlXPUWhOZgnwTFFPpCXcT\nbXHFTP7kQ9fzwfddxdO/eIMXtu6jvTN4C5CBgUF+teMgv9pxkHllxVy7dgXvunIZs0qLUlxrEUkU\nT53MZrbezPaZ2QEzu2+UfR40sxoz22lmq+OVNbPfM7M3zWzAzK4Y/6FIss2cXshHbr+ab//1h/j0\nB6877z5SACcbWvj3p1/lE5u+xxe+uZkXt+6nq7s3RbUVkUSJ26IwswzgIeB64CSwzcx+7JzbF7HP\nTcAS59wyM7sKeBhYF6fsG8D7gUcSfVAjOW9m9kR84RSWnRWgau0Kqtau4EjdWX72q928tP0g3T19\n4X12HzzJ7oMnefSHL7Fm1UKuXl3JmosXkJuTlcKai8iF8NL1tBaocc4dAzCzJ4ANwL6IfTYAjwM4\n57aaWbGZzQYWj1bWObc/tE2/tyexxRUz+fgfXsvd7/8dXt11lOptB3h9X224o6+3r5/f7DzEb3Ye\nIiuQyZqLF3D15Uu4ctVChYbIJOElKMqB2oj1OoLhEW+fco9lZQrIyc7inVcu451XLqOxpYOXth/k\nxVf3U3uqMbxPX/8Ar+w6wiu7jpAVyGT1yvmsWbWANasWUlpckMLai0gsyRrMTutWghoxyVVaXMCG\nd1/GbdddSm19E7/eeYhXdh6mtr4pvE9f/wDb3jzKtjePAsGWyZpVC7ly1QKWLpil/0YiacRLUJwA\nFkSsV4S2Re8zf4R9sj2UjWvjxo3h5aqqKqqqqsb6EZICZsaCuaUsmFvKnTe9ndr6plA31GGOR7Q0\nAI7UneVI3Vl+9LPtTCvMY/XKCi5dXsEly8t1g0IRD6qrq6murk7KZ1u8B9eYWSawn+CA9CngVeAu\n59zeiH1uBj7lnLvFzNYBX3fOrfNY9kXgfzvnto/y/S4RD9f5wjc3s/vgSQA2furW854rLRPv5Jlm\ntu8+zvY9x9h98FR4fsZI5pUVc+mKYGi8bVk5hfk5E1hTkcnJzHDOJaRpHrdF4ZwbMLN7gecIXk77\nmHNur5ndE3zbPeqce8bMbjazg0AHcHessqGDuB34JjAT+KmZ7XTO3ZSIg5L0N29WCfNmlXDrdZfS\n0dXD6/vr2L77ODv2HKe1veu8fU82tHCyoYUtL+/GgMXzy7i4ci4rK+dw0ZI5lBTlj/wlIpIQcVsU\nqaYWhb845zh0vIFdB07wxoET7D18ir7+gZhl5pYVB0Ojcg4XVc5lblmxxjjE9ya0RTFVnP+EuxRW\nRGIyM5YunMXShbO444bL6e3rZ/+R0+zaX8fr++s4XNvwljn2pxpaONXQwotb9wNQVJDL0gVlLF04\ni2ULZrF0wSzdj0pkHHwTFDI5ZWcFuGR5OZcsL+eDt15FR1cP+w7Xs+9wPXsOn6Lm2JnwU/qGtHV0\n89reWl7bO3xldtn0omBwLJzF0gVlLC6fSX5e9kQfjsikpKCQSaUgL4c1qxayZtVCIDih73DtWfYc\nOsW+w/XsPXyKzhFuG9LQ1EZDUxu/2XkovG1WaRGLK2aycN4MFlfMZFH5DMqmF6rbSiSKL4NCvwim\njuysACsr54Sf0Oec42RDCwePneHg8TPUHDvDkRPn6B9hnONMYxtnGtvYuutIeFt+bjaLymeEX/Pn\nlFI+u4SCPF1pJf7ly6CQqcvMKJ9VQvmsEq59+3IA+vsHOH6qkZpjZzh4vIGDx89Qd7p5xEtyO7t7\nw8/ciFRaXEDF7OlUzClh/pxSKuZMp2J2CdMKNfYhU59vgiLNL+6SJAoEMqmcX0bl/DJuDG3r6xug\n7nQTR0+cC75OnuVI3Tk6unpG/IzGlg4aWzrOe4gTwLTCPCpmBy/1nVtWHH7NmTmN7Czf/O8lU5x+\nksWXsrIyWVwxk8UVM8PbnHOca+7g6Mlz4QA5cbqJE2ea3zJgPqS1vYs97V1vaYEYUFpSEBEeJeEA\nUYjIZKOfVpEQM2Pm9EJmTi/kytBgOQQfznT6XCt1p5uprW+krr6J2vom6uqbRp3j4YBzzR2ca+7g\nzZqTb3l/+rR8Zs2YRllpIbNLg//OmjGNsumFlE0vIisrM1mHKTJmvgmKyCfcaTBbxiIzMyM8k3zt\nJYvC251zNDS1U3uqMTyXo/5sK6caWjhzrjXmMxWbWjtpau1k/5G3vmfA9OICykqLmBV6lZUWMqOk\nkBklBcwoKaQgL1s/xzJhfBMUIolmZuFf5NH6+wc43dgWDJAzLdSfbQmHSUNjW8wQcQyPiew/Uj/i\nPtlZAWaWFFBaUkBpcQEzS4JBUlpSwMxQmEwrzFWYSEL4JigiB7P1v44kWyCQGb76ilXnv9ffP8C5\nlg7OnGujobGN043BfxtCl+uea2qP+4T33r7+8D2wRpOZmUHptAJKpuUxfVo+00PLJUX5lEzLZ/q0\nvOC2ojwCAXV1yeh8ExQi6SIQyGT2jGnMnjFtxPdHCpKzTe00Nndwrrmds80d9PT2jVg20sDAYHii\nYTyF+TmhAMmjZFp+KGDyKS7MZVphHtMKcplWlEdxYS452Xoyod8oKETSTLwgcc7R2d0bHixvbGkP\nL59rHl4e7VLfkbR39tDe2UPd6aa4+2ZnBSguzGNaYW7olUdxYR5FBbkUF+VSVJA3HDCFueTnajxl\nsvNlUOiHViYzM6MgL4eCvBwWzC0ddb/unj4aWzpobuuiqbWT5tCrqS30b2sXza2dtLR1xu3qitTb\n1++5pQLBLrDC/ByK8nMpyM+hKD+HwoJcCvNyKCwIbi/MDy4HtwXXNWCfPnwZFCJ+kJuTFb5aK5bB\nwUFaO7qXBrY4AAAH/klEQVRpbu2ksaWTllCwNLV20NLeTVt7Ny3tXbS2d9Ha0T3qnJLRDAwM0tLW\nRUtbV/ydIxhQkJ8TDJH8XIoKcoLroZDMy82iIC+H/Lxs8nOzKcjLJj8vGDD5udnkZAcUNAnim6BI\n9+duiKRKRkZGcHyiKJ9FcR7TMtTt1dreTWt7VzBIOrpoaeumrSMiUNq7aQ0FTG9f/wXVyzHcJQat\nF3Rc+aEwycvNpjA/GCD5eTnk52aFQ6UgL5u83OzgfjlZ5OZkkZeTRV5uNnk5WWRmZlxQ/acS3wSF\niIxfZLfX3LJiT2V6+/pp7+yhraOH9s7u0C//bto7e2nv6KZtaHnovY4e2rt66BrhLsBjMTg4GBE0\nFy4QyAwGx1CI5GaTn5tFbk52KFAi3svJJi83Mmzeul9GxuQLHgWFiCRVdlaA0uIApcUFYyrX3z9A\nR1cv7V09oUAJ/tvR1UtHVw9d3X10dPXQ0dVLZ1cvnd29dHb10NEdXI/3ZMSx1KOtf4C2ju6EfF4g\nkEludoCc7AC52Vnk5GSF13OygyGTkxUgN2fk9ZzsALlDZcLvBdeTdZmzb4JCT7gTmVwCgUyKi/Iu\n+OmEfX0DdHYHQ6WzqzccIJ3dPXR29YVbLcGg6aGzu5eu7j66e/roCr26u3vHNNDvRX//AO39A+Nu\n6YwkIyMjHDqJ5JugEBF/ycrKpDjrwoMGgn9g9vT2B4Oju5ee3v5goPT00d3dR1dPb+i9oYDppbP7\n/Pe6Q+939fR5mv8yHoODg8GW1Ti77aIpKERERmFmwW6dnCymT8sf9+cNBU9Pbz/dvX109/TT09sX\nXu/p6aenb2h78L3u87YNrffTHQqeYNl+enr6GEzSRTu+DApdMiciqRAZPMUk9qFXzjn6+weDgdPb\nz5MPfiJhn+3LoBARmWrMjKysTLKyMika23UDcU2+67QukKZRiIhcGN8EhYiIXBgFhYiIxOSboNA8\nChGRC+OboBARkQvjKSjMbL2Z7TOzA2Z23yj7PGhmNWa208xWxytrZtPN7Dkz229mPzMzbzeOSQDT\nM+5ERDyLGxRmlgE8BNxI8KGOd5nZyqh9bgKWOOeWAfcAD3soez/wvHNuBfAC8PmEHNEUVl1dneoq\npA2di2E6F8N0LpLDS4tiLVDjnDvmnOsDngA2RO2zAXgcwDm3FSg2s9lxym4A/jW0/K/A7eM6Eh/Q\n/wTDdC6G6VwM07lIDi8T7sqB2oj1OoIBEG+f8jhlZzvnTgM45+rNbNZoFfjiI896qGZssR5CLyIi\no0vWzOwLGQQYdUrc9j3HxlGVt8rI0BiFiIhnzrmYL2AdsCVi/X7gvqh9Hgb+MGJ9HzA7VllgL8FW\nBcAcYO8o3+/00ksvvfQa+yve73evLy8tim3AUjNbCJwC7gTuitpnM/Ap4Admtg5ods6dNrOzMcpu\nBj4K/APwEeDHI325c05//ouIpFDcoHDODZjZvcBzBAe/H3PO7TWze4Jvu0edc8+Y2c1mdhDoAO6O\nVTb00f8A/IeZfQw4BvxBwo9ORETGzSJnLIuIiERL25nZXib5TSVmVmFmL5jZbjN7w8z+Z2j7qBMT\nzezzoUmOe83svamrfXKYWYaZ7TCzzaF1X54LMys2sx+Gjm23mV3l43PxGTN708x2mdn3zCzbL+fC\nzB4zs9Nmtiti25iP3cyuCJ2/A2b2dU9fnqjBjkS+CAbYQWAhkAXsBFamul5JPuY5wOrQciGwH1hJ\nsIvuc6Ht9wFfCi1fDLxGsPtwUeh8WaqPI8Hn5DPAvwGbQ+u+PBfAvwB3h5YDQLEfzwUwDzgMZIfW\nf0BwfNMX5wJ4B7Aa2BWxbczHDmwF3h5afga4Md53p2uLwsskvynFOVfvnNsZWm4neFVYBaNPTLwN\neMI51++cOwrU8Nb5LZOWmVUANwP/FLHZd+fCzKYB73TOfQcgdIwt+PBchGQCBWYWAPKAE/jkXDjn\nXgaaojaP6djNbA5Q5JzbFtrvcTxMdk7XoBhtAp8vmNkign85vELUxERgaGJi9Dk6wdQ6R18D/ozg\nZX5D/HguFgNnzew7oW64R80sHx+eC+fcSeArwHGCx9XinHseH56LCLPGeOzlBH+fDvH0uzVdg8K3\nzKwQ+BHwJ6GWRfTVBlP+6gMzuwU4HWphxbo8esqfC4JdB1cA/+icu4LgVYX348+fixKCf0EvJNgN\nVWBmH8SH5yKGpBx7ugbFCWBBxHpFaNuUFmpO/wj4rnNuaF7J6dB9swg1G8+Etp8A5kcUn0rn6Brg\nNjM7DPw78G4z+y5Q78NzUQfUOud+G1r/T4LB4cefi/cAh51zjc65AeAp4Hfw57kYMtZjv6Bzkq5B\nEZ7kZ2bZBCfqbU5xnSbCPwN7nHPfiNg2NDERzp+YuBm4M3TVx2JgKfDqRFU0mZxzf+6cW+CcqyT4\n3/4F59yHgJ/gv3NxGqg1s+WhTdcDu/HhzwXBLqd1ZpZrZkbwXOzBX+fCOL+VPaZjD3VPtZjZ2tA5\n/DCjTHY+T6pH8mOM8K8neOVPDXB/quszAcd7DTBA8Aqv14AdoXNQCjwfOhfPASURZT5P8GqGvcB7\nU30MSTov1zJ81ZMvzwVwGcE/nnYCTxK86smv5+KB0HHtIjh4m+WXcwF8HzgJ9BAMzbuB6WM9dmAN\n8Ebod+s3vHy3JtyJiEhM6dr1JCIiaUJBISIiMSkoREQkJgWFiIjEpKAQEZGYFBQiIhKTgkJERGJS\nUIiISEz/H9A74l46ILnTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7faa380fd090>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "thinkplot.Pdf(suite)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's how we can compute the posterior mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "333.41989326371095"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def Mean(suite):\n",
    "    total = 0\n",
    "    for hypo, prob in suite.Items():\n",
    "        total += hypo * prob\n",
    "    return total\n",
    "\n",
    "Mean(suite)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we can just use the method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "333.41989326371095"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "suite.Mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sensitivity to the prior\n",
    "\n",
    "Here's a function that solves the train problem for different priors and data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MakePosterior(high, dataset, constructor=Train):\n",
    "    \"\"\"Solves the train problem.\n",
    "    \n",
    "    high: int maximum number of trains\n",
    "    dataset: sequence of observed train numbers\n",
    "    constructor: function used to construct the Train object\n",
    "    \n",
    "    returns: Train object representing the posterior suite\n",
    "    \"\"\"\n",
    "    hypos = range(1, high+1)\n",
    "    suite = constructor(hypos)\n",
    "\n",
    "    for data in dataset:\n",
    "        suite.Update(data)\n",
    "\n",
    "    return suite"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's run it with the same dataset and several uniform priors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "500 151.849587959\n",
      "1000 164.305586423\n",
      "2000 171.338181092\n"
     ]
    }
   ],
   "source": [
    "dataset = [30, 60, 90]\n",
    "\n",
    "for high in [500, 1000, 2000]:\n",
    "    suite = MakePosterior(high, dataset)\n",
    "    print(high, suite.Mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are quite sensitive to the prior, even with several observations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Power law prior\n",
    "\n",
    "Now let's try it with a power law prior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Train2(Train):\n",
    "\n",
    "    def __init__(self, hypos, alpha=1.0):\n",
    "        Pmf.__init__(self)\n",
    "        for hypo in hypos:\n",
    "            self[hypo] = hypo**(-alpha)\n",
    "        self.Normalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what a power law prior looks like, compared to a uniform prior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGjZJREFUeJzt3W2QXFd95/Hvv580M3qWLEuWZAnFFsImwdhUbIGz6yEm\nWLAJpnarUnZ2w4ZNvN5aDK4ltWWTqpTFu6VSyQbCsti1Xggku2ZDKrFIGVAImSVeyiCM/IhkyzZ6\nlseSLI0ep6e7739f3Ds9t0czmh5P3+7T079PVVff5z59NdU/nXPuudfcHRER6U25ThdAREQ6RyEg\nItLDFAIiIj1MISAi0sMUAiIiPUwhICLSw5oKATPbZmZ7zexlM3tgivW/ZWbPJq8nzexdze4rIiKd\nYzONEzCzHPAycDtwFNgF3OXue1PbbAX2uPuImW0Dtrv71mb2FRGRzmmmJnAzsM/dD7h7BXgMuDO9\ngbs/5e4jyexTwLpm9xURkc5pJgTWAYdS84eZ+JGfyu8B336L+4qISBsVWnkwM3s/8HHgV1p5XBER\nyUYzIXAE2JCaX58sa5B0Bj8CbHP3U7PZN9lfNzESEZkld7e57N9Mc9Au4Foz22hmJeAuYEd6AzPb\nAPw18Nvu/ups9k1zd73ceeihhzpehhBeOg86FzoXl3+1wow1AXevmdl9wE7i0HjU3feY2b3xan8E\n+ENgBfAlMzOg4u43T7dvS0ouIiJz1lSfgLt/B9gyadnDqel7gHua3VdERMKgEcMBGhwc7HQRgqDz\nMEHnYoLORWvNOFisXczMQymLiEg3MDN8jh3DLb1EdK6+tfsoP3ntFJVaxEffs46br1nR6SKJiMxr\nQYXAyIUKR05dBOBcudrh0oiIzH9B9QkU8xPFqVSjDpZERKQ3BBUCpcJEcaqR+gdERLIWVAgU8hP9\nG5WaagIiIlkLKgQamoMUAiIimQssBCZqAtWamoNERLIWWAioY1hEpJ2CCoF0n8CYmoNERDIXVAik\nawJqDhIRyV6wIaCOYRGR7AUWAulLRFUTEBHJWlAhUFBNQESkrYIKgaIGi4mItFVQIdBw2wg1B4mI\nZC6oECjkNE5ARKSdggqBYkHjBERE2imsENA4ARGRtgo2BNQxLCKSvcBCQDeQExFpp6BCIJ8zLMmB\nauREerCMiEimggoBM6OQS40ViNQkJCKSpaBCAKBYSF8mqpqAiEiWgguBUsMVQqoJiIhkKbgQKOgm\nciIibRNcCOgyURGR9lEIiIj0sABDQM1BIiLtEl4IFFQTEBFpl+BCoGGcgEJARCRTwYWAxgmIiLRP\neCGQ0zgBEZF2CS8E9EwBEZG2CS8E9EwBEZG2CToE1DEsIpKtAENA4wRERNoluBAoqCYgItI2AYaA\nni4mItIuwYVASTUBEZG2CS4E1DEsItI+wYWAnicgItI+wYVA420jVBMQEclScCHQ8HjJSDUBEZEs\nBRcCDc1BqgmIiGQquBBQx7CISPs0FQJmts3M9prZy2b2wBTrt5jZD81s1Mw+PWndfjN71sx2m9mP\nZ/osjRgWEWmfwkwbmFkO+CJwO3AU2GVmj7v73tRmJ4FPAh+d4hARMOjup5opkGoCIiLt00xN4GZg\nn7sfcPcK8BhwZ3oDdz/h7k8D1Sn2tyY/B1AIiIi0UzM/zuuAQ6n5w8myZjnw92a2y8zumWljjRMQ\nEWmfGZuDWuBWdz9mZquIw2CPuz851Ybbt2+nXKnx0+eHuWrLe1h8wy1tKJ6ISHcYGhpiaGiopcc0\n98v/b9vMtgLb3X1bMv8g4O7+uSm2fQg46+5/Ms2xpl1vZu7uXChX+dTXnwGgv5Tnzz5242y/k4hI\nTzAz3N1m3nJ6zTQH7QKuNbONZlYC7gJ2XK5cqQIOmNmiZHoh8EHghct9mEYMi4i0z4zNQe5eM7P7\ngJ3EofGou+8xs3vj1f6Ima0GfgIsBiIzux+4HlgF/I2ZefJZf+nuOy9boFzqVtKR4+6YzSnoRERk\nGk31Cbj7d4Atk5Y9nJoeBq6eYtdzwLtnUyAzo5i3eqdwpeaUCgoBEZEsBDdiGHSZqIhIu4QZAuoX\nEBFpiyBDIN0vUNGdREVEMhNkCKgmICLSHkGGgJ4pICLSHkGGgJ4pICLSHkGGQPrqoDFdHSQikplA\nQyA1YEw3kRMRyUygIaBxAiIi7RB8CKgmICKSnSBDIN0xrD4BEZHsBB8CujpIRCQ7QYaAxgmIiLRH\nkCGgEcMiIu0RZAik7x2kPgERkewEGQLpmoCuDhIRyU6YIaBxAiIibRFoCKSuDlIIiIhkJtAQUHOQ\niEg7BB8CqgmIiGQn0BBINwepJiAikpUgQ6CgmoCISFsEGQJF3TZCRKQtggyBUkG3jRARaYcgQ6Ch\nOUg1ARGRzAQZAkXdSlpEpC0CDQGNExARaYfgQ0BXB4mIZCfQENCD5kVE2iHIEEh3DKtPQEQkO0GG\nwOSagLtqAyIiWQgyBMys8TnDahISEclEkCEAk54zrCYhEZFMBBsCqgmIiGQv2BDQZaIiItkLNgQK\nerqYiEjmgg2BUkNNQM1BIiJZCDYEigV1DIuIZC3YECjkUjeR051ERUQyEWwIFNUcJCKSuXBDQM1B\nIiKZCzYE0s1BqgmIiGQj2BBI1wR0iaiISDbCDQGNExARyVzAIaCOYRGRrHVJCKgmICKShWBDoOG2\nERonICKSiaZCwMy2mdleM3vZzB6YYv0WM/uhmY2a2adns+90Gm4lHak5SEQkCzOGgJnlgC8CdwDv\nBO42s3dM2uwk8Engj97CvlNquDpINQERkUw0UxO4Gdjn7gfcvQI8BtyZ3sDdT7j700B1tvtOp+G2\nEeoTEBHJRDMhsA44lJo/nCxrxlvet3HEsJqDRESyUOh0AdK2b99en169+SZgLaCrg0REAIaGhhga\nGmrpMZsJgSPAhtT8+mRZM2a1bzoEfrr/FE9/71VANQEREYDBwUEGBwfr85/97GfnfMxmmoN2Adea\n2UYzKwF3ATsus72lpme7b53GCYiIZG/GmoC718zsPmAncWg86u57zOzeeLU/YmargZ8Ai4HIzO4H\nrnf3c1Pt21TB9KB5EZHMNdUn4O7fAbZMWvZwanoYuLrZfZuxINUxfHGsNtvdRUSkCcGOGF46UKxP\nn7lY6WBJRETmr2BDYHFfYwi4q0lIRKTVgg2BUiHHQCkPQORwdnTyODQREZmrYEMAYImahEREMhV0\nCCztT4eAagIiIq3WNSEwckE1ARGRVgs7BNQcJCKSqaBDYHH/xDCGEYWAiEjLBR0CjX0CCgERkVYL\nOwQG1CcgIpKlsEMg3TGsmoCISMsFHQJLdHWQiEimgg6BxX0FLLmZ6PlyjapuKS0i0lJBh0AuZyzu\nm7hCSAPGRERaK+gQAPULiIhkKfgQ0P2DRESyE3wIaKyAiEh2wg+BAd1ETkQkK8GHgC4TFRHJTvAh\noI5hEZHsBB8CS9I3kVNNQESkpbogBNQxLCKSleBDoOEmcgoBEZGWCj4EBkp5Crn43hHlSkS5Uutw\niURE5o/gQ8DMGgaMqTYgItI6wYcA6IHzIiJZ6Y4Q0MNlREQy0RUhkL5MVFcIiYi0TpeEgPoERESy\n0BUh0NAnoOYgEZGW6YoQUE1ARCQbXRECGjAmIpKN7ggB3UlURCQTXREC6auDzl6s4u4dLI2IyPzR\nFSGwoJinrxgXtRo558u6dYSISCt0RQgArFhUqk+/cWa0gyUREZk/uiYErl4xUJ8+ePJCB0siIjJ/\ndE8IrFQIiIi0WveEwIr++vThNy92sCQiIvNH94RAqiZw6OQFokhXCImIzFXXhMCS/iLLkkFjlZrz\nxplyh0skItL9uiYEANanmoQOvql+ARGRueqqENgwqUlIRETmpqtCYL1CQESkpboqBDakxgoc0hVC\nIiJz1lUhcOWSBZQKcZFHLlT0lDERkTnqqhDI5ayhc1hNQiIic9NUCJjZNjPba2Yvm9kD02zzBTPb\nZ2bPmNmNqeX7zexZM9ttZj+ea4E1clhEpHUKM21gZjngi8DtwFFgl5k97u57U9t8CLjG3Teb2S3A\nfwe2JqsjYNDdT7WiwFc31ATULyAiMhfN1ARuBva5+wF3rwCPAXdO2uZO4GsA7v4jYKmZrU7WWZOf\n05SGy0Q1VkBEZE6a+XFeBxxKzR9Oll1umyOpbRz4ezPbZWb3vNWC1j9oeT9m8fTrI6OMVaO5HlJE\npGfN2BzUAre6+zEzW0UcBnvc/cmpNty+fXt9enBwkMHBwUu2WVDMc+WSBQyPlHGHI6cusmnVwoyK\nLiISjqGhIYaGhlp6TJvpUY1mthXY7u7bkvkHAXf3z6W2+TLwj+7+jWR+L3Cbuw9POtZDwFl3/5Mp\nPsebfWzkw99/lV2vxV0M/+bWDQxed2VT+4mIzCdmhrvbXI7RTHPQLuBaM9toZiXgLmDHpG12AB9L\nCrUVOO3uw2Y2YGaLkuULgQ8CL8ylwEDD//xfOHxmrocTEelZM4aAu9eA+4CdwIvAY+6+x8zuNbN/\nn2zzBPBzM3sFeBj4j8nuq4EnzWw38BTwLXffOddC37BhWX36xcMjlKt65rCIyFsxY3NQu8ymOQjg\nD7/5AsdOx88a/sQHruHGty3PqmgiIkFqV3NQkG7cOFEb2H3gdAdLIiLSvbo2BN6dCoFnD56mpieN\niYjMWteGwKZVC+tPGjtfrvHK8LkOl0hEpPt0bQiYGTekagPPqElIRGTWujYEYHK/wClC6eQWEekW\nXR0C77hqMf2lPAAnzo5x5JRuKCciMhtdHQKFfI5fXL+kPq+rhEREZqerQwDgxo0T4wOeeuWkmoRE\nRGah60PgXVcvpa8Yf43hkTLPHhzpcIlERLpH14dAXynfcAO5J549ptqAiEiTuj4EAD7wzisp5OKR\n06+9cZ6XX9eYARGRZsyLEFi2sMT7Nq+sz3/72WMdLI2ISPeYFyEA8MFfWlN/4tgLh89wWI+eFBGZ\n0bwJgTXL+rgpdSfRbz/7egdLIyLSHeZNCABse9ea+vSu197k4EnVBkRELmdehcCmVQu5fl08eCxy\n+OoP9lOt6UH0IiLTmVchAPBb791AIR93Dhw8eYGdzw/PsIeISO+adyGwZlkfH7lpbX1+x+6jHDut\newqJiExl3oUAwB2/tIaNVwwAUK05X/3BfiI9dEZE5BLzMgTyOePj//xt5JMBZK++cZ4dPz3a4VKJ\niIRnXoYAwPoVA3z4homrhf7umWP800vHO1giEZHwzNsQAPj1G9c23Gr6L/7fQX525EwHSyQiEpZ5\nHQL5nPEffvUarl7RD0Atcr70vVc4pPEDIiLAPA8BiO8y+qk7NrN8YfxQ+tFKxB898ZIeTC8iQg+E\nAMDyhSXuv2Nz/VGUF8o1/viJl9i9/1SHSyYi0lkWyr33zcyzLsvPj5/nC9/dx9nRavKZcPd7N/D+\n61Zh43efExHpEmaGu8/px6unQgDgjTOj/Ol39vHGmXJ92Xs2Ledjv7KRhQsKmX++iEirKATeojMX\nK3xh5z72H5/oIF6+sMjv3raJd6xdcpk9RUTCoRCYg3K1xjeeOsQP9p5IlQFu3XwF//KX17Gkv9i2\nsoiIvBUKgRbYvf8UX/2n/Zwv1+rLBkp5fuOmtbz/ulUU8j3Rdy4iXUgh0CKnz4/xtScP8NyhkYbl\nKxeV+NANa7j17VdQVBiISGAUAi32/KERHnvqIMMj5YblywaK/NovrubWt1/Boj51HotIGBQCGajW\nIv7hxTf49nOvcy65lHRcMW/88i+sYPC6VWxatVCXlYpIRykEMlSu1Pi/e4/z3eeHGblQuWT96qUL\nuOWaldxyzQpWL+3rQAlFpNcpBNpgrBrx1CsnGdpzfNpnFq9f0c+7Nyzjho3LeNsVA6ohiEhbKATa\nyN35+fHzDO05ztP7T1GuTP3s4iX9Ba5ft4Tr1savFYtKbS6piPQKhUCHlKs1njs4wo9efZMXDo1Q\nvcxTy1YuKrF5zSI2r17MNasXctWy/vrDbkRE5kIhEIDRsRovHjnDMwdP89zB0w3jDaZSKuTYsHKA\nTasWsn5FPxuvGGDN0j6NRxCRWVMIBMbdOXDiAj87eoY9R87wyvA5KrWZv1MhZ6xe2sfa5X2sXd7P\nmqV9XLWsjyuX9FEqKBxEZGqtCIHgLnr/V/d/udNFaBkHyhQYpcioFSlToGr5WR2h4BFFavHLaxSI\nKCTzORw1LInIXAQXAvOJAX1U6aMKfhGAqucoU6BsBcYozBAMRtXyVMlzcfyADWudQj0Y4nAoeEQ+\nmc8nLwWFiExHIdBm8Y/1GAt9rL4scmOMPGMUGLM8FfJUKFCxHJf88qc4RsUKNIxiuGRzJ+/eEAp5\nIvIekSMij9eX5XDVLkR6jPoEAlapRZw8O8bwmVGOnylz/GyZk+fGOHG2zImzZUanuUx1LnIGAwsK\nLFyQZ9GCQn16oFRgYEGegVKe/tLEfH8xz8CCPH3FeLnusSTSPuoY7mHuzsWxGm+eH+PkuTFOnR/j\n9PkKb54f4/SFCiMXKpy+MDbj1UqtVsgZfaU4HPqKORYk733FOCgWJNMLijkWFOL1Cwo5SoVL58eX\nlQo5DcATmYJCQGZUqUWcvVhl5GKFMxcrnB2txu8Xq5wdrXButMq5cpWzF+P36QbBdVoxb/VgKOVz\nFJPpYj7XsK6QS96TZePri/lkWT5HIZkeX1cYf8/Fywu5iWX5XPxSCEmIFALScpVaxPnRKhfGapwv\nV7lQrnGuHM9fSN4vjtXq86OVKJmPp2uXGTjXrcwgn4tDI36fCIfxZYVkvnDJ/MS28fIcuRwUcrmG\n5fWX2ZTrx5eZUT92zlLrbWK7+BiXLlOYzT9tCwEz2wb8KZADHnX3z02xzReADwHngd9x92ea3TfZ\nTiHQ5dydSs0ZrcRBUa7GATFaqVGuRJSrNUYrUcN8/B5RrtSS94ixWjw/Vo0Yq0ZNjbWQ5phBzoyc\nQa4eOvH8eJjkkiCqh4hZPQinms7VAyZ1bEvPTxwvvXyq41yyb1I2aNy/4dgW/ximj2GTPqe+H9QD\n0pLzYal96vtCw2ell6WPHx+zcT+b/PmpZa3/92zDOAEzywFfBG4HjgK7zOxxd9+b2uZDwDXuvtnM\nbgG+DGxtZl+51NDQEIODg50uxqyZGaVC3AzTisdzjp8Hd48DoRZRqaano/r0WDWiWovXVaM4WKo1\np1JLlifbV2tOJUqmI2/YplqLqETxe7XmyfqIECo3x156mqu2vGfOx3GHmjs1gC4N11adi05IBwOp\nUIqDYlKIpLYf3xcag60VmrlE9GZgn7sfiAtijwF3Aukf8juBrwG4+4/MbKmZrQY2NbGvTNKtIdBq\n4+fBzOIO4+JsBtq1ThRNBMJ4cFSjOChqk9ZNbOvxj22ybRRBNYrq249vV5vmFfnEMaLIOTz0Ate9\nfzBZR8N2teR4tfFpj8ucXu/OZe9x1U26OQTGQ7hBh8O4mRBYBxxKzR8mDoaZtlnX5L4iQcvljFLO\nOnoLj5M/XMnvf3jLnI/jSSjU3PGIieAYDw+PgyNK1kVJkIwHS33aJ+2TmvYkpNLHcuL1PulYXv9M\nGubHg8vdceLjOeCRc2LlAO/bvDKeT44fJduNH2P8czx1bKBezvHjul+6T1RfN7584ryQ7DPe9xWl\njjN+zChKbZfaJ1RZDRZT75NIgMzizupuHiX62j8s49/dtqnTxZi1dDDVgwOHVKikgwniQBlfBhOh\nQ7Lvo/fMvVwzdgyb2VZgu7tvS+YfjL/PRAevmX0Z+Ed3/0Yyvxe4jbg56LL7po4RcFaKiISpHTeQ\n2wVca2YbgWPAXcDdk7bZAXwC+EYSGqfdfdjMTjSxLzD3LyIiIrM3Ywi4e83M7gN2MnGZ5x4zuzde\n7Y+4+xNm9mEze4X4EtGPX27fzL6NiIjMSjCDxUREpP06frcvM9tmZnvN7GUze6DT5WknM1tvZt83\nsxfN7Hkz+1SyfLmZ7TSzl8zsu2a2tNNlbRczy5nZT81sRzLfk+ciucz6r8xsT/L3cUsPn4v/ZGYv\nmNlzZvaXZlbqlXNhZo+a2bCZPZdaNu13N7PPmNm+5O/mg818RkdDIDWY7A7gncDdZvaOTpapzarA\np939ncB7gU8k3/9B4HvuvgX4PvCZDpax3e4Hfpaa79Vz8XngCXe/DriBeGxNz50LM1sLfBK4yd3f\nRdyEfTe9cy6+Qvz7mDbldzez64HfBK4jvnvDl6yJYcqdrgnUB6K5ewUYH0zWE9z99fHba7j7OWAP\nsJ74HPx5stmfAx/tTAnby8zWAx8G/kdqcc+dCzNbAvwzd/8KgLtX3X2EHjwXiTyw0MwKQD9whB45\nF+7+JHBq0uLpvvtHgMeSv5f9wD6aGJfV6RCYbpBZzzGztwHvBp4CVrv7MMRBAVzZuZK11X8F/jPx\n5dLjevFcbAJOmNlXkqaxR8xsgB48F+5+FPhj4CDxj/+Iu3+PHjwXKVdO890n/54eoYnf006HgABm\ntgj4JnB/UiOY3Fs/73vvzexfAMNJzehyVdh5fy6ImzxuAv6bu99EfMXdg/Tm38Uy4v/5bgTWEtcI\n/jU9eC4uY07fvdMhcATYkJpfnyzrGUkV95vA19398WTxcHLvJcxsDfBGp8rXRrcCHzGz14D/Dfyq\nmX0deL0Hz8Vh4JC7/ySZ/2viUOjFv4sPAK+5+5vuXgP+BngfvXkuxk333Y8AV6e2a+r3tNMhUB+I\nZmYl4sFkOzpcpnb7n8DP3P3zqWU7gN9Jpv8t8PjkneYbd/8Dd9/g7r9A/HfwfXf/beBb9N65GAYO\nmdnbk0W3Ay/Sg38XxM1AW82sL+nkvJ34woFeOhdGY+14uu++A7gruXpqE3At8OMZD97pcQIWP2/g\n80wMJvsvHS1QG5nZrcAPgOdJbhsC/AHxP9z/IU71A8BvuvvpTpWz3czsNuD33f0jZraCHjwXZnYD\ncQd5EXiNeABmnt48Fw8R/8egAuwGfg9YTA+cCzP7X8AgsBIYBh4C/hb4K6b47mb2GeB3ic/V/e6+\nc8bP6HQIiIhI53S6OUhERDpIISAi0sMUAiIiPUwhICLSwxQCIiI9TCEgItLDFAIiIj1MISAi0sP+\nP/3fX0VG0qG+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7faa0a9f9b50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "high = 100\n",
    "hypos = range(1, high+1)\n",
    "suite1 = Train(hypos)\n",
    "suite2 = Train2(hypos)\n",
    "thinkplot.Pdf(suite1)\n",
    "thinkplot.Pdf(suite2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's see what the posteriors look like after observing one train."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEPCAYAAACQmrmQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW9+P/Xe2ay7wkQlkBAFqm0VawiV/1q1N4qbrS2\nvQVvr5V6W6rV2uX2Z6v1Frra1rbqtVVbqdW6gLd2wWoVrQavOyqIVZYgECBAWEIgZJ+Z9++PcxJm\nJstMICeTTN7Px2MeOcvnc877EzTvOZ9zzucjqooxxhjjBV+yAzDGGJO6LMkYY4zxjCUZY4wxnrEk\nY4wxxjOWZIwxxnjGkowxxhjPeJ5kROQCEVkvIhtF5IYeytwhIlUiskZEZkZsXyIitSKyNqb8iSLy\nioisFpHXReQUr9thjDGm7zxNMiLiA+4EzgdmAPNFZHpMmTnAZFWdCiwE7orYfZ9bN9ZPge+q6kzg\nu8DPPAjfGGPMMfL6SmYWUKWq1araDiwF5saUmQs8AKCqrwEFIlLqrr8IHOjmuGGgwF0uBGo8iN0Y\nY8wxCnh8/HHA9oj1HTiJp7cyNe622l6O+zXgaRH5OSDA6cceqjHGmP42VG/8Xw1cr6oTcBLO75Ic\njzHGmG54fSVTA0yIWC+ja9dWDTA+TplYn1PV6wFU9Y8isqS7QiJiA7MZY8xRUFXpj+N4fSWzCpgi\nIuUikg7MA5bHlFkOXAEgIrOBelWN7CoT9xOpRkTOduucB2zsKQBVTdnPd7/73aTHYO2ztln7Uu/T\nnzy9klHVkIhcC6zASWhLVHWdiCx0dutvVPVJEblQRDYBjcCCjvoi8jBQAZSIyDacJ8ruA74I3C4i\nfqDFXTfGGDPIeN1dhqo+BRwfs+2emPVre6h7eQ/bXwI8fTdmd30LL2zYy+RRuXxkUpGXpzLGmJTl\neZIZqu5duZmte5uAWhZfNoNxxVnJDqmLioqKZIfgqVRuXyq3Dax95gjp7/63wURE9Gjb95/3vtG5\nfOVZEzlz2oj+CssYYwY1EUH76ca/XckYYzwxceJEqqurkx2G6UV5eTlbt2719ByWZIwxnqiuru73\nJ5VM/xLpl4uVXg3VlzGNMcYMAZZkjDHGeMaSTAJCIbvkN8aYo2FJJgFvbu1uIGhjzHB09dVX88Mf\n/rBz/a677mL06NHk5+dz4ID9rYhljzD3IPIR5lMmFfGl8yb3V1jGDAvuY7DJDsNTwWCQ/Px8Xn/9\ndT74wQ8mO5w+6+nfqD8fYbYrGWOMOUq7d++mtbWVD3zgA0dVP9WTMFiSMcYMQz6fj82bN3euL1iw\ngP/+7/8GYOXKlYwfP55f/OIXlJaWMm7cOH7/+993KVtVVcX06c5Ev0VFRXz0ox8F4OWXX2bWrFkU\nFRVx2mmn8corr3TWPeecc/jOd77DmWeeSU5ODlu2bOGcc87h5ptv5owzziAvL4+5c+dSV1fHZz/7\nWQoKCjjttNPYtm3bAPxWvGHvyRhjBtwnr7+7X4/32O1f6lP5eO+H7N69m4aGBnbu3MmKFSv41Kc+\nxSc+8QkKCgo6y0ydOpV3332X4447joMHDyIiHDhwgIsvvpg777yTefPm8eijj3LRRRfx/vvvU1Tk\njIH44IMP8tRTTzFt2jTC4TAAy5YtY8WKFZSUlDB79mxOP/107rrrLh544AEWLFjA4sWLWbKk2xlN\nBj27kjHGDDvxuqnS09O5+eab8fv9zJkzh9zcXDZs2BD3eE888QTTpk3j8ssvx+fzMW/ePKZPn87j\njz/eWfbKK69k+vTp+Hw+AgHne/6CBQuYOHEieXl5zJkzh8mTJ3POOefg8/n49Kc/zerVq/uh1clh\nSSYBm/c2JjsEY8wAKikpwec78ucxOzubw4cPx623c+dOysvLo7aVl5dTU3NkHsbx48fHVqO0tLRz\nOSsrq8t6IucerKy7LAF1h9uSHYIxKaWv3Vv9LTs7m6amps713bt3d/vHv6/Gjh3LY489FrVt27Zt\nzJkzp3N9IIZyGUzsSiYB2en+ZIdgjOlHM2fO5OGHHyYcDvPUU0+xcuXKoz5WZNfbhRdeSFVVFUuX\nLiUUCrFs2TLWrVvHJZdc0h9hD0meJxkRuUBE1ovIRhG5oYcyd4hIlYisEZGZEduXiEitiKztps51\nIrJORN4RkVu8bIMxJrXcdtttLF++nKKiIh555BE+8YlP9Fq+t6uPyH3FxcX87W9/49Zbb2XEiBHc\neuutPPHEE503/bs7Tqpf2Xj6MqaI+ICNwHnATmAVME9V10eUmQNcq6oXichpwO2qOtvddyZwGHhA\nVT8cUacCuBG4UFWDIjJCVfd1c/5+eRkzO93PHVfM7KW0MSbWcHgZc6hLhZcxZwFVqlqtqu3AUmBu\nTJm5wAMAqvoaUCAipe76i0B34zRcDdyiqkG3XJcEY4wxJvm8TjLjgO0R6zvcbb2VqemmTKxpwFki\n8qqIPC8ipxxzpL1obg95eXhjjElZQ/XpsgBQpKqzReRU4FHguO4KLlq0qHO5oqLiqObmVoVQWPH7\nUrvv1BgzPFVWVlJZWenJsb1OMjXAhIj1MndbbJnxccrE2g78CUBVV4lIWERKVHV/bMHIJNMXPoFw\nRFdl7cEWxhZlHdWxjDFmMIv9Ar548eJ+O7bX3WWrgCkiUi4i6cA8YHlMmeXAFQAiMhuoV9XaiP3i\nfiL9BTjXrTMNSOsuwRyLtED0r8ZuXxpjTN95mmRUNQRcC6wA3gWWquo6EVkoIl90yzwJbBGRTcA9\nwDUd9UXkYeBlYJqIbBORBe6u+4DjROQd4GHcJGWMMWZw8fyejKo+BRwfs+2emPVre6h7eQ/b24H/\n6K8YjTHGeMPe+DfGGOMZSzIJ2neoNdkhGGOGkUmTJvHcc88lO4xjZkkmQTX1zckOwRhjhhxLMgkK\n2DsyxhiPpPLwO5ZkjDHDzqRJk7jllluYMWMGJSUlXHXVVbS1HZnS47e//S1Tp05lxIgRfPzjH2f3\n7t2A897dV77yFQCCwSC5ubnccIMz7m9LSwtZWVnU19cD8Oqrr3LGGWdQVFTEzJkzo0Z67m4a5t6s\nWrWK008/naKiIsaNG8d1111HMBjsc0zJMFTf+DfGDGGRA9D2h3v/s+8jSz388MM888wzZGdnc/HF\nF/ODH/yA733vezz33HPceOONPPvss5xwwgl84xvf4DOf+QwrV67k7LPP5qtf/Srg/OEfPXo0L7zw\nAgAvv/wy06dPp7CwkJqaGi6++GIeeughzj//fP7xj3/wyU9+kg0bNlBSUgJET8Mc70rG7/dz2223\nceqpp7J9+3bmzJnDr3/9a77yla8kHFOy2JVMglL4ataYYem6665j7NixFBYWctNNN/HII48ATvK5\n6qqrOPHEE0lLS+PHP/4xr7zyCtu2beNf/uVfqKqq4sCBA7zwwgtcddVV1NTU0NTUxAsvvMDZZ58N\nwEMPPcRFF13E+eefD8B5553HKaecwpNPPtl5/shpmP3+3uesOvnkk5k1axYiwoQJE/jiF7/YeWWU\naEzJYkkmQSvX7012CMaYflRWVta5XF5ezs6dO4GuUyjn5ORQUlJCTU0NmZmZnHLKKVRWVvLCCy9Q\nUVHB6aefzosvvth5pQNQXV3No48+SnFxMcXFxRQVFfHSSy91drtB99Mw96SqqopLLrmEMWPGdCbF\nffucwecTjSlZrLssQSPyMpIdgjEp42i6t/rb9u1HBn+vrq5m7NixgDOFcnV1dee+xsZG9u/fz7hx\nzuDwZ511Fs899xxr1qzh1FNP5ayzzuLpp59m1apVnHXWWYCTQK644gruuSfqvfMofZms7Oqrr+bk\nk09m2bJlZGdnc/vtt0dN85xITMliVzIJsmfLjEktv/rVr6ipqaGuro4f/ehHzJs3D4D58+dz3333\nsXbtWlpbW7nxxhuZPXs2EyY4Y/2effbZPPDAA5xwwgkEAgEqKiq49957mTRpUuf9ls9+9rM8/vjj\nrFixgnA4TEtLCytXruy8WuqrhoYG8vPzyc7OZv369dx1111R+xOJKVksyRhjhqXLL7+cj33sY0yZ\nMoWpU6dy0003Ac79k+9///tcdtlljBs3ji1btrB06dLOeqeffjotLS2d3VAnnHACWVlZUd1SZWVl\n/PWvf+VHP/oRI0eOpLy8nFtvvZVwOAwkdhUTWebWW2/loYceIj8/n4ULF3YmxL7ElCyeTr+cbMcy\n/fKX73+L1vZw5/qMcfl8bc60/grNmJQ3mKdfnjRpEkuWLOHcc89NdihJlQrTL6eMd2sOJTsEY4wZ\ncizJ9MFg/VZmjOmbvtx0N8fGni4zxgw7mzdvTnYIw4ZdyRhjjPGM50lGRC4QkfUislFEbuihzB0i\nUiUia0RkZsT2JSJSKyJre6j3DREJi0ixV/EbY4w5ep4mGRHxAXcC5wMzgPkiMj2mzBxgsqpOBRYC\nkQ+A3+fW7e7YZcC/AtXd7TfGGJN8Xt+TmQVUqWo1gIgsBeYC6yPKzAUeAFDV10SkQERKVbVWVV8U\nkfIuR3X8EvgmsNy78KM1tATJz0obqNMZM6SVl5fbDfZBLnL4HK94nWTGAdsj1nfgJJ7eytS422p7\nOqiIXApsV9V3BvI/4up9TXxofMGAnc+YoWzr1q3JDsEMAkPu6TIRyQJuxOkq69zcU/lFixZ1LldU\nVFBRUeFVaMYYMyRVVlZSWVnpybG9TjI1wISI9TJ3W2yZ8XHKRJoMTATeFucypgx4U0Rmqeqe2MKR\nScYYY0xXsV/AFy9e3G/H9vrpslXAFBEpF5F0YB5d76EsB64AEJHZQL2qRnaVCRFXKqr6T1UdrarH\nqeoknC64md0lGGOMMcnlaZJR1RBwLbACeBdYqqrrRGShiHzRLfMksEVENgH3ANd01BeRh4GXgWki\nsk1EFnR3GgZokOR3dxwciNMYY0zK8PyejKo+BRwfs+2emPVre6h7eQLHP+6YAuwDG1TGGGP6xt74\n7wN7GNMYY/rGkowxxhjPWJLpg5aI+WWMMcbEZ0mmD17cuC/ZIRhjzJBiSaYPSgsykh2CMcYMKZZk\n+sDGYTLGmL6xJGOMMcYzlmT6YHd9S7JDMMaYIcWSTB81t4WSHYIxxgwZlmT6qLE1mOwQjDFmyLAk\nk4DsdH+yQzDGmCHJkkwi7KEyY4w5KpZk+mhfQ2uyQzDGmCHDkkwCmlqP3OzffdCeMDPGmERZkknA\n+OKszmWxvjNjjEmYJZkEjCnMil/IGGNMF54nGRG5QETWi8hGEbmhhzJ3iEiViKwRkZkR25eISK2I\nrI0p/1MRWeeWf0xE8r1uR4f39xweqFMZY8yQ52mSEREfcCdwPjADmC8i02PKzAEmq+pUYCFwV8Tu\n+9y6sVYAM1T1JKAK+LYH4XdqajvybkxdY5uXpzLGmJTi9ZXMLKBKVatVtR1YCsyNKTMXeABAVV8D\nCkSk1F1/ETgQe1BVfVZVOyZ3eRUo8yh+AMaXZHcu2zszxhiTOK+TzDhge8T6Dndbb2VquinTm88D\nfz+q6BI0uiDTy8MbY0zKCiQ7gGMhIjcB7ar6cE9lFi1a1LlcUVFBRUXFMZ3zra31x1TfGGMGm8rK\nSiorKz05ttdJpgaYELFe5m6LLTM+TpkuRORK4ELg3N7KRSaZo+WzeWSMMSks9gv44sWL++3YXneX\nrQKmiEi5iKQD84DlMWWWA1cAiMhsoF5VayP2CzEDu4jIBcA3gUtV1fNX8KeOzu1cTg/YU9/GGJMo\nT/9iqmoIuBbnabB3gaWquk5EForIF90yTwJbRGQTcA9wTUd9EXkYeBmYJiLbRGSBu+t/gFzgGRF5\nS0R+7WU70vyWWIwx5mh4fk9GVZ8Cjo/Zdk/M+rU91L28h+1T+y3APmoLhlFVm4rZGGMSYF/RE+CL\n+S3VHGhOTiDGGDPEWJJJQFpMlrHZMY0xJjGWZBI0pTQ3fiFjjDFRLMkchZ0HbLh/Y4xJhCWZBEWO\nWVZ7yJKMMcYkwpJMgiaOODJ+mT3SbIwxibG/lgmaEDFIpt34N8aYxFiSSZBGLL+wYW/S4jDGmKEk\noSQjIn8SkYvc+WGGJb/vyMuXZcU2U6YxxiQi0aTxa+ByoEpEbhGR4+NVSDXTx+QlOwRjjBlyEkoy\n7iRh/w6cDGwFnhWRl0VkgYikeRngYBE5jMzWvU1JjMQYY4aOhLu/RKQEuBL4T2A1cDtO0nnGk8gG\nmcy06F9Va7vd/DfGmHgSGiBTRP6MM8jlH4BLVHWXu2uZiLzhVXCDyZjC6Psw7SElY1hcwxljzNFL\ndBTm37pD8ncSkQxVbVXVUzyIa1DKyfDT2GpXMMYYk6hEu8t+0M22V/ozkKFAI55j3r7f7ssYY0w8\nvSYZERktIh8BskRkpoic7H4qgOze6kYc4wIRWS8iG0Xkhh7K3CEiVSKyRkRmRmxfIiK1IrI2pnyR\niKwQkQ0i8rSIFCQSy7FqingJ81Bz+0Cc0hhjhrR4VzLnA7cCZcAvgJ+7n68DN8Y7uPtezZ3ucWYA\n80VkekyZOcBkdyKyhcBdEbvvc+vG+hbwrKoeDzwHfDteLP3hI5OKOpe1l3LGGGMcvd6TUdX7gftF\n5JOq+thRHH8WUKWq1QAishSYC6yPKDMXeMA932siUiAipapaq6ovikh5N8edC5ztLt8PVOIkHk/5\nIx5jfv39OmZPKfH6lMYYM6T1mmRE5LOq+iAwUUS+HrtfVX8R5/jjgO0R6ztwEk9vZWrcbbW9HHeU\nqta6MewWkVFx4ugX9U1HRmLOzvAPxCmNMWZIi9ddluP+zAXyuvkMFgPSexXZXWaMMSa+eN1l97g/\nFx/l8WuACRHrZe622DLj45SJVdvRpSYio4E9PRVctGhR53JFRQUVFRXxo+5BZtqRq5dXN9XxnxXH\nHfWxjDFmsKisrKSystKTY8frLrujt/2q+pU4x18FTHHvq+wC5gHzY8osB76M82LnbKC+oyusIwz3\nE1vnSuAnwOeAv/YUQGSSOVZZ6UeSTI51lxljUkTsF/DFi4/2uqKreC9jvnksB1fVkIhcC6zA6Zpb\noqrrRGShs1t/o6pPisiFIrIJaAQWdNQXkYeBCqBERLYB31XV+3CSy6Mi8nmgGvi3Y4kzUZNH5Ua0\nbSDOaIwxQ1siT5cdE1V9CmdImsht98SsX9tD3ct72F4HfPRYY+uriIfLaGoLEQ4rPl/sRZYxxpgO\n8brLblPVr4rI43Rzc11VL/UsskEo8p4MwM76ZsqKE3on1RhjhqV43WV/cH/e6nUgQ0F6IPphvLB1\nmRljTK/idZe96f5cKSLpwHScK5oNqtrWW91UNb44i+11zQBU72tkQoldyRhjTE8SHer/IuBu4H2c\nJ70michCVf27l8ENRrvqWzqXDzUHkxiJMcYMfomOwvxz4BxVrVDVs4FzgF96F9bgdcpxR17IbG6z\nYf+NMaY3iSaZBlXdFLG+GWjwIJ5BrzA7vXP51U37kxiJMcYMfvGeLrvMXXxDRJ4EHsW5J/NpnBct\nh7XSgsxkh2CMMYNavHsyl0Qs13Jk5OO9QFbX4qlvRlk+T63dDcCGXcPyYs4YYxIW7+myBb3tH47S\n/DGPMdsLmcYY06NEny7LBK7CmXiss49IVT/vUVyD1nEjc6LWQ6r4ugytZowxBhK/8f8HYDTOLJUr\ncUZKHpZ9RT6fEIi4crEnzIwxpmeJJpkpqnoz0OiOZ3YRcJp3YQ1uwYhX/bfsaUxiJMYYM7glmmTa\n3Z/1IvJBoAAYkNkoB6PIYf4bWu2FTGOM6UmiSeY3IlIE3Iwzl8t7OMPtD0sTI+7LvLfjUBIjMcaY\nwS2hG/+qeq+7uBIY9tNBtgXDnctZNnmZMcb0KKErGREpEZH/EZG3RORNEblNREq8Dm6wOvW44s7l\nlev2JjESY4wZ3BLtLlsK7AE+CXwK2AcsS6SiiFwgIutFZKOI3NBDmTtEpEpE1ojISfHqisiJIvKK\niKwWkddF5JQE29EvIof8L8lN76WkMcYMb4kmmTGq+n1V3eJ+fgCUxqskIj7gTpxHn2cA80VkekyZ\nOcBkVZ0KLMQZ7Tle3Z/iTMU8E/gu8LME29EvppYemYZ5/+FhOeOBMcYkJNEks0JE5omIz/38G/B0\nAvVmAVWqWq2q7ThXRHNjyswFHgBQ1deAAhEpjVM3jPOEG0AhUJNgO/pF7AyZ9Y2WaIwxpjvxBshs\nwBkQU4CvAg+6u3zAYeC/4hx/HLA9Yn0HTvKIV2ZcnLpfA54WkZ+7sZ0eJ45+VZCdFrXe3B6icCAD\nMMaYIaLXKxlVzVPVfPenT1UD7senqvkexZTIGC1XA9er6gSchPM7j2LpUXHEvZiNuw4P9OmNMWZI\nSOgRZgARuRQ4y12tVNW/JVCtBpgQsV5G166tGmB8N2XSe6n7OVW9HkBV/ygiS3oKYNGiRZ3LFRUV\nVFRUJBB2fHUR92Ia2+yFTGPM0FVZWUllZaUnxxZVjV9I5BbgVOAhd9N84A1V/Xacen5gA3AesAt4\nHZivqusiylwIfFlVLxKR2cBtqjq7h7rzVHW9iLwLXKOqK0XkPOAWVT21m/NrIu3rzpfvf4vWdud9\nmDuvmElmevR9mAdfqqbSfXx5bFEm3/vkB4/qPMYYM9iICKraLyP/JnolcyFwkqqG3QDuB1YDvSYZ\nVQ2JyLXACpyuuSWquk5EFjq79Teq+qSIXCgim4BGYEEvdde7h/4CcIebiFqAL/ahzf0i8oXMnIyE\nLwiNMWZY6ctfx0Kgzl0u6K1gJFV9Cjg+Zts9MevXJlrX3f4yMKDvxsQ6eWIRL1c50y9X7bZ7MsYY\n051Ek8yPgdUi8jzOjfmzgG95FtUQMNqmXjbGmLjiJhkREeBFYDbOfRmAG1R1t5eBDXYj8qLf9G9u\nC5GVbuOYGWNMpLhJRlVVRJ5U1Q/hjMBsAH/MlMtb9zXygbFePdVtjDFDU6Jv/L8lIl2e3hrORCRq\nXpn3amzIf2OMiZXoPZnTgM+KyFacJ8AE5yLnw14FNhSkB3w0tjrTL9s0zMYY01WiSeZ8T6MYomZN\nLubptbUAVK7by2fPKE9yRMYYM7jEG7ssE/gSMAV4B+ddFXu93VWcY8P8G2NMb+Ldk7kf532Ud4A5\nwM89j2gImTo6L2o9GAr3UNIYY4aneN1lJ7hPleGOD/a69yENHWMLo9+VaWgJUmRXN8YY0ynelUx7\nx4J1k3UV8Ef/+qpq7c1/Y4yJFC/JnCgih9xPA/DhjmURsWd2Y6yzx5iNMSZKr91lqmqvsMfx4fEF\nrN1+ELCpmI0xJlaiL2OaHkwcmdO5bC9kGmNMNEsyx2ja6Nxkh2CMMYOWJZljNKEkO2r9QKN1mRlj\nTAdLMscoO2bCso27GpIUiTHGDD6eJxkRuUBE1ovIRhG5oYcyd4hIlYisEZGTEqkrIteJyDoReced\nHjpp8jKPJJoNuy3JGGNMB0+TjIj4gDtxxj6bAcwXkekxZeYAk1V1KrAQuDteXRGpAC4BPuS+LHqr\nl+2IJ7LL7M0tB5IYiTHGDC5eX8nMAqpUtVpV24GlwNyYMnOBBwBU9TWgQERK49S9Gril4wVRVd3n\ncTt6NTXi5n/HqMzGGGO8TzLjgO0R6zvcbYmU6a3uNOAsEXlVRJ4XkVP6Neo+OnliUdR6W9DGMDPG\nGEh8qP+BJPGLEACKVHW2O5nao8Bx3RVctGhR53JFRQUVFRUJBbG3roHq3fXkZGUQ1jDQ83upY4uy\notZ31DVx3Ch7tNkYMzRUVlZSWVnpybG9TjI1wISI9TJ3W2yZ8d2USe+l7g7gTwCqukpEwiJSoqr7\nYwOITDJ98f72vYRUaG0L8srqzZw3+/iE6764cZ8lGWPMkBH7BXzx4sX9dmyvu8tWAVNEpFxE0oF5\nwPKYMsuBKwBEZDZQr6q1cer+BTjXrTMNSOsuwRwL1SPLf3p2ddzyE0ceufn/jjvMjDHGDHeeXsmo\nakhErgVW4CS0Jaq6TkQWOrv1N6r6pIhcKCKbcKZ2XtBbXffQvwN+JyLvAK24ScorO/fGTxozxhWw\ndW8TAAca2+OUNsaY4cHzezKq+hRwfMy2e2LWr020rru9HfiPfgwzrvb2EGlpPd+XOXliIU+s2dW5\n3hoMkRGw8UWNMcObvfHfA4l5/OC1tVt6LV9WHD28zPqd9lKmMcZYkumBxDzktuzvq3ot7/dFl39p\nY1Jf3THGmEHBkkyCErkvE/nm/1tb670MxxhjhgRLMn3Q0tr7Df0zjx8xQJEYY8zQYEmmD5b9/Y1e\n958xtSRq3Yb9N8YMd5Zk+mD582/3uj8j5umzp9bu9jIcY4wZ9CzJ9FEo1Pu4ZNnpRxLNP97d43U4\nxhgzqFmS6aOnXny31/3nnDBqgCIxxpjBz5JMH/3uTy/1uv+s6dE3/3fXt3gZjjHGDGqWZBIw64Pl\nUesaObBZjJLcjKj1v6/d1UNJY4xJfZZkEnDFx/8lav0v/1jTa/nIFzNf2tiv43YaY8yQYkkmAdmZ\n6VHrDz7+Wq/lLzslel623q58jDEmlVmSSdAFZ86IWm9sbu2xbMUHRkatv73Nhv43xgxPlmQS9PnL\nzohaX/Srv/VYNvZ9mYdervYkJmOMGewsySTI74/+VW3evrfX8lNKj8yMafPLGGOGK0syffDday6O\nWv/Hq+t6KAmXnz4han3PIXuU2Rgz/HieZETkAhFZLyIbReSGHsrcISJVIrJGRE5KtK6IfENEwiJS\n7GUbOnz4+LKo9V8/srLHspEjMgPc/Y/NnsRkjDGDmadJRkR8wJ3A+cAMYL6ITI8pMweYrKpTgYXA\n3YnUFZEy4F+BAb3hce5pUeFTvbPnR5TTA0d+vdv2N3kWkzHGDFZeX8nMAqpUtdqdMnkpMDemzFzg\nAQBVfQ0oEJHSBOr+Evimx/F3cc38s6PWv/6T/+2x7NXnTY5aty4zY8xw43WSGQdsj1jf4W5LpEyP\ndUXkUmC7qr7T3wHHIyKUFOZEbavZ0/0EZR8sy49a/9Uz73sWlzHGDEaBZAfQDel1p0gWcCNOV1nc\nOosWLepcrqiooKKi4tiiA+78znzm/9e9netf+eFSHrv9S93FSnrAR1vQGbm55kDzMZ/bGGP6W2Vl\nJZWVlZ5bJLlEAAAZsklEQVQc2+skUwNEPmZV5m6LLTO+mzLpPdSdDEwE3hYRcbe/KSKzVLXL2PqR\nSaa/pKcFyMpMp7nlyKRkm6r3MKW86wjM/3XhNH60fH3n+oZdDRw/Jq/fYzLGmKMV+wV88eLF/XZs\nr7vLVgFTRKRcRNKBecDymDLLgSsARGQ2UK+qtT3VVdV/qupoVT1OVSfhdKPN7C7BeOm+H3wuav2G\nX/yp23LHjcqNWv/ZExs8i8kYYwYbT5OMqoaAa4EVwLvAUlVdJyILReSLbpkngS0isgm4B7imt7rd\nnYY4XWxeSEvzM2FM9JPTPU3PHPliJkAwzsRnxhiTKiSVB28UET3a9p3y9ccIulWf/96FFOVldSmj\nqnzqq/dEbVv28y8QCEQPK3O4JchXHzwycvPJEwu55qNTjiouY4zxmoigqv3y5d3e+D8GIsKVHz89\nattnvvHbLuVyM6Nvfb21tfun0YwxJtVYkjlGl5zz4S7bXlrd9VHl6z4WfeXy2iabZ8YYk/osyfSD\n+34Y/RDAL37/DG3twahtJ04ojFr/beUWz+MyxphksyTTD/JzszjzI9FXKpHv0XQ46/gRUeubag97\nGpcxxiSbJZl+8rUrPtpl2+e/c3/U+n+cWR61fsvj6zHGmFRmSaYfPfyzq6LWDzY0s/z5tzvXRYSJ\nI6NHZ964u2FAYjPGmGSwJNOPMtLT+N51l0Ztu/8vr7B+8+7O9Rsv+UDU/p/+zV7ONMakLksy/WzG\nlLFUzDo+attNt/+FXXsPAuDzSZeXM/9vQ++zbBpjzFBlScYD1/37OV22XfuDR9hT53SN3XBxdBK6\n//+qSeWXYo0xw5clGY90Nyrz1YsfYkftAUSEj32oNGrfTf/7z4EKzRhjBowlGQ91l2iu/9Ey1qzf\nzr+dNj5q+55DrRxobOtS3hhjhjJLMh77420Lu2z7/l1P8JtH/4+bPx79EMA3H1k7UGEZY8yAsCTj\nMRHpNtE8/dK7fP37f8AZRPqI7/35vQGKzBhjvGdJZgCISLddZwB177xJKGLo/237m9i8x0YCMMak\nBksyA+ix27/ERWd/KGqbALXvvUP1ziMDZv5o+frOKZuNMWYosyQzwD5/2Rn84oZPR23Lop26ukO8\n8W41obCTXK75/Vv2WLMxZsjzPMmIyAUisl5ENorIDT2UuUNEqkRkjYicFK+uiPxURNa55R8TkXyv\n29GfyseWdLlPM1GdK5nV67bz5nvVAHxhyZsDHpsxxvQnT5OMiPiAO4HzgRnAfBGZHlNmDjBZVacC\nC4G7E6i7ApihqicBVcC3vWyHFzru01z20Zmd2yaFnTf/VeGNd6vZUXuAq++zRGOMGbq8vpKZBVSp\narWqtgNLgbkxZeYCDwCo6mtAgYiU9lZXVZ9V1Y6bFq8CZR63wzP/fslpLPv5FwDnH2NCuK5z3+59\nh3hl7Vb+446VSYrOGGOOjddJZhywPWJ9h7stkTKJ1AX4PPD3Y440iQIBP4/d/iVuvvoi0ggxOnww\nav87W/ZxylcfZcljLyYpQmOMOTqB+EUGnCRcUOQmoF1VH+6pzKJFizqXKyoqqKioOJbYPHXS9PE8\ndvuX+OUDz/LEmzXs8x0ZSDMofv7nxVqefOFuwHnJUyThX5UxxvSosrKSyspKT44tXj7BJCKzgUWq\neoG7/i1AVfUnEWXuBp5X1WXu+nrgbGBSb3VF5ErgC8C5qtraw/n1aNt3ytcfI+hWff57F1KUl3VU\nxzkWn1n0GBsOdt1+XHhvZyb+ydcvY0r5qAGNyxiT2kQEVe2Xb7Fed5etAqaISLmIpAPzgOUxZZYD\nV0BnUqpX1dre6orIBcA3gUt7SjCpYNmiT/I/XzqLwtzoBLfZN5Kgm2Zu+MWf+OT1d/Od2//KgUNN\nyQjTGGN65OmVDHQmhNtxEtoSVb1FRBbiXJX8xi1zJ3AB0AgsUNW3eqrrbq8C0oGONxhfVdVrujn3\nkL6S6VC1u4Gf/G0D++sPs6XmyEubhdpEiTZ2Kf+haeP48vwKRhbnDWSYxpgU0Z9XMp4nmWRKlSQD\n0NDSztcedKZyDofDvL99LwcPtwDR3Wex8nOz+PLlFXzkhAl2D8cYkxBLMglKpSQDoKpdXtBsammj\namstmcEmRmlD3GOcNH08n73kNCaOK7GkY4zpliWZBKVakunw4EvVVK7rOmVzQ1ML/j3VtDbETzYd\npk0s5VMfO5mTpo/H77dRhowxlmQSlqpJBqC5LcR1D6zucf/155Xz57+/xtsbdvT52OefMYPzzzyB\nCWOK7WrHmGHIkkyCUjnJdPjzGzU8sWZXj/u/fcl0Jo7I5pU1m/ndn1/iYENzn88xoiiXD08r47QT\nJ3H8xFLycjKPJWRjzCBnSSZBwyHJQPf3amKdOW0EV541EYBgMMSrb29h2d9XsXNvNy/iJKgj+Zw4\nvYzJ40cyekS+XfkYkwIsySRouCSZDqGwsvB3vSebUfkZfPWCqYzKj74aqdlTz9MvvstLb71PfcOx\nv29TVlrEB6eOZdrEUsrHFjN2VCHpaYNxgAljTCxLMgkabkmmQzisfDFOsgGYNjqXqyomUZKb0WWf\nqrKnroE3/lnN6nXbWL1uezdHODq52RmUjS5iwphixo8uYuyoQsaOKqSkIMcePjBmELAkk6DhmmQi\n3f50Fe9sj98llpcZ4HP/byInTijotcsrGAyxtWY//9y0k39W1fD2+h2EPfhvqCAvi7zsTMaPKWbC\nmCJGFuVRmJ/NqJI8ivOzycpM7/dzGmMclmQSZEnmiLrDbdz5zCa27U+sK2za6FzmfmQc00bnJnyf\n5WBDM5t37KOqupbqnXVs27n/mO75JMLv95EW8JObndF5Xyjg9zOutJCCvCzyczIpzM8mLzuTtDS/\np7EYkyosySTIkkz3tuxt5FfPbKK+qT3hOnmZAf71Q6WcPqWEwpy+X0WEw2Fq9zdQs6eebTvr2Lar\nju27D7Br70Fa2xKPoz9kZqThE8HnE06cPp66+kY+OG0sjU2tHFc2ElVlZHEemRkBcrMzycvJJCsj\nzbryzLBhSSZBlmTi23OohYde3sa7Ow71uW75iGxOKi/kA2PzKR+RTdox/hEOBkPsPXCYuoON7N53\nkO27DrC3roHDza3s2nuQ/fVdx2kbaD4RwqqkpwUoHZFPU3Mr0yaOZm9dAzOmjGFffSNTJ4yiPRii\npDCHnGznfldhbhZpaQFys9PJSE8jMz1AIGBXVmZwsiSTIEsyfaOqvLn1AH9eVUPtoaMf3Hp0YSZT\nS3OZUprL5FG5jMhLJ9CPVwGqSlNLG4cOt7C//jB76w6zp66BsCqbt+/lYEMzPp9QVb0Hn89HOByO\nf9AkyUhPo7WtnTEjCwgGw+w90MCHp5VR39BEelqACWOKaWpuJSMjjdKSfNrag2RlplOUn0UopORk\nZ5CTlU5awI+IkJOVQXqan7SAn8yMtM5le7Tc9IUlmQRZkjk2qkpV7WGe+Wctq7fW98sxi3LSGFuU\nxdjCLMYVZzG6IJMRuenkZ6Xh83n3h7AjMTU0ttLQ2Mz++kbCqmzbVYdPhOqa/YjPx669B2lra+dQ\nYwuHm1JrFom0gJ+szHRaWttpaw8yYUwxaQE/h5taaW5tZ1r5KAJ+H/vqG8nPzWRkUR5pAT9NLW34\nfMKoEmddFYKhECMKc/H7fQT8PoLBMAV5WQT8PgIB5wtFZobz5SIt4O8slxbwE/D78Pt9lvgGMUsy\nCbIk0/+CoTBb9zWxdls9q6vr2VXf0u/nyE73U1acRW5mgNyMAKUFmYzMz6AgK428rAB5mWlkpg3s\nHylVpaW1ncbmNhqb22hpbeNQYwtt7SH2HTiMT4Rtu+ooyM1kw9ZaCvOzaWhsYcuOfZSNLmL95t0U\n5Wd3zvkjQOr+n5cYv99HKORcZebnZpEW8OH3+Qirsu/AYSaVjXCSlt+P3y+89/5uTjx+HBlpAXxu\n0mptC3LwcDMTx5Z0lgv4/Rw83Ex2Zjr5uZn4/c5x/T4fre1BMtIC5GSl4/P78Im4+4W29hB5ORlO\nWXdfIOBDFdLTA1FlC3KzUrq705JMgizJDJyWthDb6prYsKuB9Tsb2L6/iaa20ICcuzA7DREoyE6j\nKCcdvwhjCjNRhTGFmYgIhTlppPl95GYGyAj4yE7392sXXl+pKsFgmObWNlrbgrQFQxxubCEUVg42\nNBMKhzsflvD7fChK9c46sjPT8ImPDVtryc/JJD3Nz+59h6jdf4iy0UW0t4fYUXuA/NwsgsEQTS1t\nSWtjKivIy+I7Cy/kuPEjkx2KJ/ozyXj+CrY78dhtHJl47CfdlLkDmIMzadmVqrqmt7oiUgQsA8qB\nrcC/qaq3z8qaXmWm+5k2Oo9po/O4ZGb0vrZgmL0Nrew80MyOumZqDjSzq76ZvYdaCffDd5yOp+QO\nNLazde/RjVZQWpBB7cFWppTmEvALOw80M6OsgLZgmLAq4wqzCIYVn0BJbgYBv9AWDFOcm975wENO\nhp80v4+AX0jz+9yPON1DMV2BIkJamp+0NO+/vKgq7cEQbe0h2oMhWlrbCYbCBIMhWtuCNLc6v7/2\nYIgDBxsJBHyEw06dvXWH8fsEf8BPKBhi9/4G2tuD5Odm0R4M0d4eZPX6HUyfVEooHKY9GGbDlt0U\n5mWTk5XunCcUYn99I2nuN/9QWAf1fbJEHGxoZuWqqpRNMv3J0yQjIj7gTuA8YCewSkT+qqrrI8rM\nASar6lQROQ24G5gdp+63gGdV9acicgPwbXfbsFJZWUlFRUWyw4grPeBjXFEW44qyOPW47su0BcPU\nN7Wxv6GNusY2Dja380LlSvImnkhTa5C0gI8ddX0f3DNRtQed+y+bag93bnul6sgspKs59ntSaW7y\n8fuEjW+/xgdmziY7w0/A5+Nwa5CG5namjs7DJxDw+2hsDdLaHqasOIuAXzoT1d5DrZSPyOks5xMI\nhpVQWCnKTsfnc56CC/icx7Rb2kPkZzr3vPw+QUTw+wL409LwB5SSAj8+cepMGC+d9Z0P+HzOst/n\nrIu7veNY3Yn336aqEgqFaWsPEQqHO5NeMBSmrT1IOOzsD4bChMJhDh1uIT3NTyisBEMhwiGlsbmV\nlrYgmemBznKhUJidew6Sm51BwO9ztymhcJjqnfvJycogKzPNTXROsjt0uIWtO/czZcJIQqEwIffc\noXCYLTv2kZmRRlF+NqFQ2O0ubWXvjo20BU845v8mhgOvr2RmAVWqWg0gIkuBucD6iDJzgQcAVPU1\nESkQkVJgUi915wJnu/XvByqxJDOkpQd8jMrPjBpT7fU/r+em6+d1KauqtIeUwy1BDrcGaWwNcqCx\njVBY2VHXTHa6nw27GhiVn8k/dxxkVH4GG3cfJifDT2Or04Xn9wmh/riM6oP2kNIecs6/a8ObjDn+\nI13eVXqvpuuj5N29QPvPo3jk3Gtpfiep+UV47c+P8C87CjuTl4jzRaKxNURpfgZIdCITt4wv4mdb\nMMy2/U1MH5uHACJBRJz7WU65dEQykJaIbT5oL87h/cY2yoty3HoQQJhSWubUF+nc7nOT5Ifcn85+\n2N/QRl5WgDPSA86ss+55126o4fnX17Nzx/uejHSRirxOMuOAyEGvduAknnhlxsWpW6qqtQCqultE\nRvUUwKd++ORRBR60/34GLREhPSAU56ZTnHv0w8sEQ2Ga20Odf/yCoTCtwTAHm9rx+4QDjW0EQ4qI\nc6XjE0gL+Khvamd3fQulBRm0h5R3th9k4shswmHYXteEKuRnBQiGdMDuSw0G7SGFkHYuH24Jdluu\nr4/Hr9+Z+CR8kXYe6P+HUvbUNVInORyWDLYdGj7/tsdiMA6LezQ3m3pMCZv2HXsXi6+HLgEztAX8\nPvLc+yklud6dR9XpygqF1b2aCfPj90fx9U99kFBYaQuFaWkPITj3K4JhpaUtRFso3HnFFQorre1h\n6pvaycsKdG4Lh5W9Da0IQnaGn3DEuVqDYTbvaWTSyByna0iVsNK5f1PtYcqKs5zzuvsiy9UdbiMz\nzeneU3XLhBVVOssMZ3XNw/wXkChV9ewDzAaeilj/FnBDTJm7gc9ErK8HSnurC6zDuZoBGA2s6+H8\nah/72Mc+9un7p7/ygNdXMquAKSJSDuwC5gHzY8osB74MLBOR2UC9qtaKyL5e6i4HrgR+AnwO+Gt3\nJ++vR/CMMcYcHU+TjKqGRORaYAVHHkNeJyILnd36G1V9UkQuFJFNOI8wL+itrnvonwCPisjngWrg\n37xshzHGmKOT0i9jGmOMSa6UHLtcRC4QkfUistF9j2bIEZEyEXlORN4VkXdE5Cvu9iIRWSEiG0Tk\naREpiKjzbRGpEpF1IvKx5EWfGBHxichbIrLcXU+ZtgG4j+P/rxvzuyJyWqq0UUS+JiL/FJG1IvKQ\niKQP5baJyBIRqRWRtRHb+tweETnZ/Z1sFJHbBrodPemhfT91418jIo+JSH7Evv5rn5c3/pPxwUmc\nm3BGA0gD1gDTkx3XUbRjNHCSu5wLbACm43QV/n/u9huAW9zlE4DVOF2gE93fgSS7HXHa+DXgQWC5\nu54ybXPj/j2wwF0OAAWp0EZgLLAZSHfXl+HcGx2ybQPOBE4C1kZs63N7gNeAU93lJ4Hzk922Xtr3\nUcDnLt8C/NiL9qXilUznC6Cq2g50vMQ5pKjqbnWH11HVwzhP1JXhtOV+t9j9wMfd5UuBpaoaVNWt\nQBVd30kaNESkDLgQuDdic0q0DcD9Vvj/VPU+ADf2g6ROG/1AjogEgCyghiHcNlV9ETgQs7lP7RGR\n0UCeqq5yyz0QUSepumufqj6rqh3j+7yK8/cF+rl9qZhkenq5c8gSkYk430JeJeZFVKDjRdTYdtcw\nuNv9S+CbOI9LdkiVtoEzYsU+EbnP7RL8jYhkkwJtVNWdwM+BbThxHlTVZ0mBtsUY1cf2jMP5e9Nh\nKP3t+TzOlQn0c/tSMcmkFBHJBf4IXO9e0cQ+qTHkntwQkYuAWvdKrbfHzIdc2yIEgJOBX6nqyThP\nTn6L1Pj3K8T5ll+O03WWIyL/Tgq0LY5Uaw8AInIT0K6qj3hx/FRMMjXAhIj1MnfbkON2RfwR+IOq\ndrwLVOuO7YZ7+brH3V4DjI+oPpjbfQZwqYhsBh4BzhWRPwC7U6BtHXYA21X1DXf9MZykkwr/fh8F\nNqtqnaqGgD8Dp5MabYvU1/YMuXaKyJU43daXR2zu1/alYpLpfAFURNJxXuJcnuSYjtbvgPdU9faI\nbR0vokL0i6jLgXnuUz6TgCnA6wMVaF+o6o2qOkFVj8P593lOVf8DeJwh3rYObjfLdhGZ5m46D3iX\nFPj3w+kmmy0imSIiOG17j6HfNiH6yrpP7XG71A6KyCz393IFPbwoniRR7RNnKpVvApeqauSAcv3b\nvmQ/9eDRkxQX4DyNVQV8K9nxHGUbzgBCOE/HrQbecttVDDzrtm8FUBhR59s4T4KsAz6W7DYk2M6z\nOfJ0Waq17UScLz1rgD/hPF2WEm0EvuvGuRbnpnjaUG4b8DDOlCKtOEl0AVDU1/YAHwHecf/23J7s\ndsVpXxXOy+xvuZ9fe9E+exnTGGOMZ1Kxu8wYY8wgYUnGGGOMZyzJGGOM8YwlGWOMMZ6xJGOMMcYz\nlmSMMcZ4xpKMSTkiEhaRn0Wsf0NE/rufjn2fiFzWH8eKc55Pich7IvKPmO3lIhI7u2yix3yxf6Iz\nJnGWZEwqagUuE5HiZAcSSUT8fSh+FfCfqnpezPZJRA8BkvDxVfXMPpzfmH5hScakoiDwG+DrsTti\nr0REpMH9ebaIVIrIX0Rkk4j8WEQuF5HXRORtd3iNDv8qIqvEmRjvIre+z50E6jV3EqgvRBz3BRH5\nK86wMrHxzHcngVorIj92t92MM//HEhH5SUyVHwNnuiM7Xy8inxORv7pXPM+KSI6IPCsib7hxX9pD\nW5+XIxOq/SGizC3iTEa2RkR+2qffujHdCCQ7AGM8oMCvgHe6+SPdXdkOH8aZGK4eZ1Ku36rqaeLM\nSnodR5JWuaqeKiJTgOdFZDLO2Fb1bvl04CURWeGWnwnMUNVtkScWkTE4k0XNdM/5jIhcqqrfF5Fz\nga+r6uqYeL8FfENVL3WP8Tm3/odU9aCI+ICPq+phESnBmR6iY+y+yLaehDM51W431tOB9W7d6e6x\n8zHmGNmVjElJ6kyLcD9wfR+qrVLVParaBryPM14VOGM1TYwo96h7jk1uuenAx4ArRGQ1zuyBxcBU\nt/zrsQnGdSrwvDqjGYeBh4CzIvb3Ng1CpGfUmRANnP+nfywib+OMuzVWREZ1U+d1Vd2lzrhSa9z2\nHQSaReReEfkE0Jzg+Y3pkSUZk8pux7m3kROxLYj73707kmx6xL7IkWjDEethoq/6I68IxF0X4DpV\nnel+JqszkRc4c8n0JNFE0pvI4/87MAKYqaozcYanz+ymTmRbQ0BAnWH7Z+FML3Ex8FQ/xGaGOUsy\nJhUJgKoewLnquCpi31bgFHd5Ls7owX31aXFMxrkRvwF4GrhGnDmAEJGp4syE2ZvXgbNEpNi9aT8f\nqIxTpwHI62V/AbBHVcMicg7OxGIdek1obryFqvoUTtfgh+PEYkxcdk/GpKLIK42fA1+O2PZb4K9u\nt9bT9HyV0dvw5NtwEkQesFBV20TkXpwup7fcK6Q9xJn/XFV3i8i3OJJY/qaqf4tz/rVA2I3/93Sd\nl/4h4HG3u+wNnKHa47WpY3s+zu+m48rna73Fb0wibKh/Y4wxnrHuMmOMMZ6xJGOMMcYzlmSMMcZ4\nxpKMMcYYz1iSMcYY4xlLMsYYYzxjScYYY4xnLMkYY4zxzP8PCHbB1h5fa7gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7faa0aa81250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dataset = [60]\n",
    "high = 1000\n",
    "\n",
    "thinkplot.PrePlot(num=2)\n",
    "\n",
    "constructors = [Train, Train2]\n",
    "labels = ['uniform', 'power law']\n",
    "\n",
    "for constructor, label in zip(constructors, labels):\n",
    "    suite = MakePosterior(high, dataset, constructor)\n",
    "    suite.label = label\n",
    "    thinkplot.Pmf(suite)\n",
    "\n",
    "thinkplot.Config(xlabel='Number of trains',\n",
    "                 ylabel='Probability')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The power law gives less prior probability to high values, which yields lower posterior means, and less sensitivity to the upper bound."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "500 130.708469863\n",
      "1000 133.275231375\n",
      "2000 133.997463081\n"
     ]
    }
   ],
   "source": [
    "dataset = [30, 60, 90]\n",
    "\n",
    "for high in [500, 1000, 2000]:\n",
    "    suite = MakePosterior(high, dataset, Train2)\n",
    "    print(high, suite.Mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Credible intervals\n",
    "\n",
    "To compute credible intervals, we can use the `Percentile` method on the posterior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(69, 869)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hypos = xrange(1, 1001)\n",
    "suite = Train(hypos)\n",
    "suite.Update(60)\n",
    "\n",
    "suite.Percentile(5), suite.Percentile(95)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have to compute more than a few percentiles, it is more efficient to compute a CDF.\n",
    "\n",
    "Also, a CDF can be a better way to visualize distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEPCAYAAABRHfM8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYFNXVx/HvQTZRUNxwAVHBJcGgoqCihsFEQSW4K+CK\ne4wagxvmjYq+JopL1MT4uuFGENwj7ogyGqMCgoiJbEpENsWoEESWAc77RxXQPTPdUwNdXb38Ps/D\nQ99bt6pOlzhnqm7de83dERERqUuDpAMQEZHioIQhIiKRKGGIiEgkShgiIhKJEoaIiESihCEiIpHE\nmjDMbIiZfWVmk7O0+ZOZzTCzSWa2d5zxiIjI+ov7DuNhoEemjWZ2BNDO3XcFzgfujTkeERFZT7Em\nDHd/B/guS5OjgcfCtmOBzcysVZwxiYjI+km6D2MHYHZKeW5YJyIiBSbphCEiIkWiYcLnnwu0SSm3\nDutqMDNNeiUish7c3XJxnHwkDAv/1GYk8CvgCTM7AFjo7l9lOpAmSgwMGjSIQYMGJR1GQdC1WEfX\nYp1SuRbuzpSZX/LUqxOYPH1OvfZt2HAjenX7CacffWDO4ok1YZjZ40AFsKWZfQFcBzQG3N3vd/eX\nzexIM/sUWAL0jzMeEZFC5u7869N5PPHKB3zy2fx67bvl5ptw3M870a3zrmzctPHa+tNzGF+sCcPd\n+0Voc1GcMYiIFCp3Z+rMLxnxynj+OWNevfZts21Ljj+sEwfuvQsNG24UU4Tpku7DkPVQUVGRdAgF\nQ9diHV2LdQr1Wrg7M2f/hyde+YAJn8yq177td9yG4w/vxH4ddqRBg2TeV7Ji6RcwMy+WWEVE1pi3\nYCFPvjqBv0+YUa/9dtupFSf22Jd9ftQGs/Xvszazour0FhEpG4sWL+W50R/yQmXGGZFqtWvbbTix\nx750+vGOG5Qg4qSEISKyAZavqOL1d6fw5KsTWLJ0eeT92m6/JSf22Jf9O+6U2COm+lLCEBGpB3dn\n4idf8PhL4/l87n8i77fNFs05qed+HLJv+7x1UueaEoaISB3mf72I4S+P5x8TP428z8ZNG3Nij33p\ncdCPadqkUYzR5Y8ShohINSuqVvLy2/9k6Mj367XfkT/dk+MP78TmzZvFFFmylDBEpOy5O598Np/H\nnn+fT79YEHm/fX7UhlN/sT9tt9+yYDuqc0kJQ0TK0vc/LOe50R/ytzcmRd5n65bNOeOYA4uqozqX\nlDBEpGx8PH0uDz37D76Y/23kfU44vBO/6L4XmzZrEmNkxUEJQ0RK1rLlVTwzaiLPjv4w8j4dd2vN\nqb/Yn3Y7bh1jZMVJCUNESsq/5/yHh579R+TJ+5o0bkS/ozrT46AONGpUnK+75osShogUtdWrV/Pm\n2Gnc/9TfWbVqdaR99vlRG844pitttm0Zc3SlRQlDRIrO4iXLGPbiWF5/d0rkffof21V3ERtICUNE\nisKsed/wwFPvMGVmtEdNe+yyLf2P6Ur7ttvEHFn5UMIQkYLk7oz/5yzuGV7J4iXLIu3Tq1tHTj5i\nP5pt3LjuxlJvShgiUjBWrlzFK3//F4/87d1I7Zs2acQvT+7GQZ3alcXAuaQpYYhIopYtr+Lp1ybw\nXMQBdHvssi3nHH8wO7feKubIpDolDBHJu8VLlvHwc+/y1vjpkdpXdNmdM44+gBabbhxzZJKNEoaI\n5MU3C7/n4efe471Jn0Vq3+fIzhxz6N56q6mAKGGISGy+Wfg99wx/i0lTZ9fZtkGDBlzYpxvdOu9a\nlvM0FQMlDBHJqXkLFnLx70dEatt8k6ZcfEr3gl6WVNZRwhCRDbbg28X88vphkdpuu1ULLj7lUPbY\nZduYo5JcU8IQkfWyaPFSzvrdo5Hatm7VkotP6a5BdEVOCUNEIlu2vIpTrhwSuf0NF/emQ/vtY4xI\n8kkJQ0SycncGDH4q8hoSv7vgKPb5UZuYo5IkKGGISK3uf/LvvPaPf0Vqe3n/wzlw711ijkiSpoQh\nImu9M+FT7nhsdKS25514CD0O7hBzRFJIlDBEytz8rxdx0Y3DI7Xt1a0j/Y/rGnNEUqiUMETK0MqV\nqzj96kdYvqKqzra777wtv//10RonIUoYIuXk8RfH8czrEyO1HXbL2TRt0ijmiKSYKGGIlLhZ875l\nwOAnI7W9Y+BJ7LjdFjFHJMVKCUOkBLk7Z1/zGIsWL62zrTqvJSolDJES8vYH07lr6Jt1tmvXZmsG\nX3ac+iWkXpQwRIrciqqV9L38wUht77/+VLbcfNOYI5JSpYQhUqReGDM50lKmZx9/EEf+9Cd5iEhK\nXewJw8x6AncCDYAh7j642vYWwF+BHYGNgNvd/ZG44xIpRstXVNHvirrncmrQoAHDbz2bhg21+JDk\njrl7fAc3awBMB34GzAPGA33cfWpKm6uBFu5+tZltBUwDWrn7ymrH8jhjFSlkY8ZO4+7Hx9TZ7toL\ne7HX7q3zEJEUCzPD3XPSWRX3HUYXYIa7zwIwsxHA0cDUlDYONA8/Nwe+qZ4sRMrR6tWrOXnAA6yu\n4xel7bfejD/9Tx91YEvs4k4YOwCpazPOIUgiqe4GRprZPGBT4OSYYxIpaJ/OWsBVf3y2znb/e8nR\n/LjddnmISCRQCJ3ePYAP3f1QM2sHvG5mHd39++oNBw0atPZzRUUFFRUVeQtSJG5/HjaGynHTsrZp\nsenGPHTj6bqbkIwqKyuprKyM5dhx92EcAAxy955heSDgqR3fZvYicJO7/yMsvwFc5e4fVDuW+jCk\n5ER9Jfay/ofRde92eYhISk0x9WGMB9qbWVtgPtAH6FutzSzg58A/zKwVsBswM+a4RBIV9bGT5nOS\nQhJrwnD3VWZ2ETCKda/VTjGz84PNfj9wI/CImU0Od7vS3aMt7SVSZIa/NI6nR2Wf/O/nB/6IX/bp\nlqeIRKKL9ZFULumRlBQrd+eca4aycPEPWdsNHnAc7dtuk6eopFwU0yMpkbK1dNkKTr3qoTrbjbjt\nXBo10gA7KXxKGCI5NnfBQi75/Yisbbr8ZCeuOqdnniISyQ0lDJEcGffx5wx+8NWsbQaceRgH7aO3\nnaQ4KWGIbKCnR01k+Evjsra559p+tNqyRZ4iEomHEobIerpr6Bu8/cGMrG2G33YOjRvpfzMpDfqX\nLFJPA//4LDNmLci4fdNmTXjkD2dqNLaUHCUMkQjcnXOvHcp3/838aqw6sqXU1ZkwzGwC8BDwuLt/\nF39IIoXD3el7+YNUrVyVsc0JPfal75Gd8xiVSDLqHLhnZu2B/gSzyH4APAyMyvcoOg3ck3xyd04a\n8ACrV6/O2Oaift3pvv/ueYxKpP5yOXAv8kjvcDGkXsD/AasIEsdd+ZrGQwlD8sHdOeXKh1i+oipj\nm4Hn9qTznjvlLyiRDZD3kd5m1pHgLuNI4BlgGHAw8Cawdy4CEUnaBYOG8fV3izNuH/SrX/CT3XbI\nY0QihSVqH8ZCYAgw0N2Xh5vGmtlBcQYnkg+/u+t5psycn3H7DRf3pkP77fMYkUhhitKHsYu7z6xW\nt7O7/zvWyGrGoUdSklN/ebySN8dOzbhddxRSCvLah2FmE929U7W6Ce6+by4CiEoJQ3LlxcrJPPzc\nuxm3X33eEezXoW0eIxKJT176MMxsD6ADsJmZHZeyqQXQNBcnF8mnj6bN4YZ7Xsy4/ZJTD6Vb593y\nGJFIccnWh7E7wVtRmwO/SKlfDJwbZ1AiufT1t4u54PphGbf3PaoLJxzeKeN2EQlEeSR1oLu/l6d4\nssWhR1JSL3Wtl911n3ZcduZheYxIJP/y0odhZle6+y1m9megRiN3vyQXAUSlhCFRuTsnXnpfzX+0\noc2bN2PIjafnNSaRpORrHMaU8O8PcnEikXx4+Nl3efGtyRm3P3XHeTRo0CCPEYmUDq3pLSVh2r+/\n5Ld3/i3j9kdv6s+mzZrkMSKRwpCvt6ReoJZHUWu4e+9cBCCyIZavqKLfFUMybr/pN8ey206t8hiR\nSOnK9kjqtrxFIbIezvrdoyxavLTWbScfsR8n9dwvzxGJlDY9kpKi89b46fzpr2/Wuq1J40Y8fuvZ\neY5IpHDl65HUk+5+kpl9TPqjKQPc3TvmIgCRqJYsXc7pAx/OuH3YLWfTtEmjPEYkUl6yvVa7nbvP\nN7Na50hw91mxRlYzHt1hlLHzrhvKNwuX1Lrt+ot+wZ67as4nkdrkfT0MM9sW6EJwpzHe3b/Mxcnr\nQwmjPE2eNofrM0znsdfurbn2wl55jkikuOR78sFzgGsJ1r4woBtwg7s/lIsAolLCKC+rVq3mpAH3\nZ9w+4rZzadRoozxGJFKc8p0wpgFd3f2bsLwl8K6753VtSiWM8jHshbE8O/rDWrdpynGR+sn3invf\nEEw4uMbisE4kp7J1arfasgX3XNsvzxGJSKpsb0kNCD9+SrC63vMEfRhHA5nnXhBZD7+/72UmfvJF\nrdse+cOZNN9EM+qLJC3bHUbz8O/Pwj9rPB9fOFJuvvvvD5xzzWO1bjvh8E70PapLniMSkUw0cE8S\n8z93/Y2pM2t/4e6J28+lYUN1aotsqLz2YZjZ1sCVBKvvrX0u4O6H5iIAKT+LFi/lrN89Wuu2S0/7\nGYfst2ueIxKRKKJ0eg8DniBYfe8C4Azg6ziDktJ160OjeP+jmbVu09TjIoUtymu1E9x9XzObvGY6\nEDMb7+6d8xLhujj0SKqILVtexSlX1j6r7BVnHc4Be+2S54hEykO+X6utCv+eb2ZHAfOALaKewMx6\nAncCDYAh7j64ljYVwB1AI+Brd+8e9fhS+F4YM5lH/vZurdt0VyFSPKLcYfQC/g60Af4MtACud/eR\ndR7crAEwHfgZQaIZD/Rx96kpbTYD3gUOd/e5ZraVu/+nlmPpDqPIuDsnXHpfrdt+1beCQw/YI88R\niZSfvN5huPuaiXwWAfX9zb8LMGPNRIVmNoJgHMfUlDb9gGfcfW54vhrJQorPzNlfc8Vtz9S6TdN6\niBSnKG9J7QLcBRwIrAbeA37j7rX3XKbbAZidUp5DkERS7QY0MrMxwKbAn9x9aIRjS4G67Jan+Xxu\nzbzf46AOnHfSIQlEJCK5EKUP43HgL8CxYbkPMBzYP4cxdAIOBTYB3jOz99z90xwdX/JkRdVK+l7+\nYK3bhtx4Ops3b5bniEQkl6IkjGbVfuP/q5ldEfH4c4EdU8qtw7pUc4D/uPsyYJmZvQ3sRTAlSZpB\ngwat/VxRUUFFRUXEMCRuH06ZzY33vlTrtqfvPB+znDxCFZE6VFZWUllZGcuxsy2gtOZNqKuA74AR\nBHNJnQy0dPer6zy42UbANIJO7/nAOKCvu09JabMHQWd6T6AJMBY42d0/qXYsdXoXqAsGDePr7xbX\nqB94bk8677lT/gMSkbXy1ek9gSBBrDnR+SnbHKgzYbj7KjO7CBjFutdqp5jZ+cFmv9/dp5rZawQT\nGq4C7q+eLKQwVVWtos/lD9S6bfht59C4UZQbWBEpFppLStZLpregttt6M+7+Xd8EIhKR2uR7LqlG\nwC+Bn4ZVlcB97l6VcScpaf834i1GvzelRv21F/Zir91bJxCRiORDlIF7DxKMwF4zW9xpwCp3Pyfm\n2KrHoTuMhGUbiPf4rWfTpHGjPEckInXJ9xKtH7n7XnXVxU0JI1nZ5oLSW1AihSvfc0mtMrN27v5Z\nePJdCDqnpUx89sXXXHl7zf6KM4/pyi+6d0wgIhFJQpSEcQUwxsxmErwx1RboH2tUUjCGvzSOp0dN\nrFH/l2v6se1WLRKISESSkjVhhJMHLgV2BXYPq6e5+/K4A5Pk9bnsAapW1ryZ1Gp4IuUpSh/Gh+6+\nT57iyRaH+jDyZPXq1Zz4m/tr3ab+CpHikss+jCgLEbxhZsebfkqUhWXLq2pNFr2778Uzd12gZCFS\nxqLcYSwmmBRwJbCMoB/D3T2vD7B1hxG/Bd8u5pfXD6tRf8PFvenQfvsEIhKRDZXv9TCa5+JEUtg+\nmjaHG+55sUa9ZpkVkTUyJgwz2wb4LdCeYJ6nm939v/kKTPLnb29MYujI92vUq3NbRFJlm632VYIJ\nCN8GegHN3f3M/IVWIx49korBzQ+8yvh/fl6jXp3bIqUhLyO9q4/mNrOJ7t4pFyddH0oYuXfaVQ/x\nw7IVaXWNGm7EiNvPTSgiEcm1vPVhmFlL1k1vvlFq2d2/zUUAkozjf31vjbpD9t2VS0//WQLRiEgx\nyHaH8TnBGt61ZSZ3911ijKu2eHSHkQOZJhDsf2xXelVomg+RUpPXyQcLhRLGhsuULH573hHs26Ft\nAhGJSNzyPfmglIBMo7dvu+IEdm69VQIRiUixUcIoAytXruLky2oupXrPtf1otaUmEBSRaJQwSlym\nZKEBeSJSX1HmksLMDjaz/uHnrc1s53jDklzIlCwevam/koWI1FuUuaSuA/YDdnf33cxse+Apdz8o\nHwGmxKFO73rIlCyG3nwWzTZunEBEIpKEfM9WeyzQG1gC4O7zAM0vVcBWrVpda7IYdsvZShYist6i\nJIwV4a/2DmBmm8QbkmwId+ekATXfhhp2y9k0bdIogYhEpFRESRhPmtl9wOZmdi4wGqj566skLtM4\ni78OPkvJQkQ2WKSBe2Z2GHA4wajv19z99bgDqyUG9WHUobbpPh67uT+bbNwkgWhEpBDkdaS3mQ0A\nnnD3ubk44fpSwsiutmTx8O/PoMWmGycQjYgUinx3ejcHRpnZ383sIjNrlYsTS+6cUEuyuG/QqUoW\nIpJTdSYMd7/e3TsAvwK2A94ys9GxRyaR3HDPi1S/77rz6pPZquWmicQjIqUr0sC90ALgS+AbYJt4\nwpH6GPHKeD6aNiet7oaLe9Nm25YJRSQipazOhGFmF5pZJfAGsCVwrrtrHuyEjf/n5zz16oS0uktO\nPZQO7bdPKCIRKXVR5pJqA1zq7pPiDkaimbtgITc/8GpaXe/ue9Gt824JRSQi5SDbAkot3P2/ZrZF\nbdvzveKe3pIKLFm6nNMHPpxWt0ubrbn18uMTikhEClm+1vR+0d17mdm/CUZ5p55QK+4lINPAvGfu\nuiCBaESkGGjFvTJV21gLJQsRySav4zDM7I0odRKv864bWqPuyT+el0AkIlKuMnZ6m1lToBmwlZm1\nZN0jqRbADnmITULPjf6QbxYuSat7+PdnsNFG9XkrWkRkw2T7iXM+MAHYI/x7zZ/ngbujnsDMeprZ\nVDObbmZXZWnX2cyqzOy4qMcuB7PmfcNfXxibVveHS4/RKG4Rybsoc0ld7O5/Xq+DmzUApgM/A+YB\n44E+7j61lnavA0uBh9z92VqOVXZ9GMtXVNHviiFpdScc3om+R3VJKCIRKTa57MOocxyGu//ZzPYE\nfgw0Tal/LMLxuwAz3H0WgJmNAI4GplZrdzHwNNA5YtxloXqyMFCyEJHE1JkwwiVaKwgSxsvAEcA7\nQJSEsQMwO6U8hyCJpB5/e+AYd+9uZvppGDrtqodq1D2tN6JEJEFRek1PIHik9KW79wf2AjbLYQx3\nAql9Gzm5dSpmo9+bwg/LVqTVDb/tnISiEREJRJkaZKm7rzazlWbWgmASwjYRjz8X2DGl3DqsS7Uf\nMMLMDNgKOMLMqtx9ZPWDDRo0aO3niooKKioqIoZRPL5dtIT/G/FWWt3tV55A40ZR/lOJSLmrrKyk\nsrIylmNH6fS+B/gt0Ae4DPgemBTebdS170bANII7lPnAOKCvu0/J0P5h4IVy7fSubSR3r24d6X9c\n14QiEpFil+9O7wvDj/ea2atAC3efHOXg7r7KzC4CRhE8/hri7lPM7Pxgs99ffZd6xF5yapv2Q8lC\nRApFtrmkOmXb0d0nxhJRBqV+h/Hm+1P5y/DKtLqn7zyf4EmdiMj6ydcdxu1ZtjlwaC4CEPj+h+U1\nksU91/ZTshCRgpIxYbh793wGUs7OuDp9uvLe3fei1ZYtEopGRKR2UcZhnF5bfcSBe1KHG+99qUbd\nGcccmEAkIiLZRXlXM3X0dVOCN54mEm3gnmQx+8vv+HDK7LS6Ebedm1A0IiLZRXlL6uLUspltDoyI\nLaIy4e5cetMTaXXXXdiLRo02SigiEZHs1md+7CXAzrkOpNycfU36DVqzpo3puHvrhKIREalblD6M\nF1g3PqIBwZxST8YZVKmb/vlXLFq8NK1u6OCzEopGRCSaKH0Yt6V8XgnMcvc5McVT8tydq+94Lq3u\nrt+enFA0IiLRRenDeAsgnEeqYfh5C3f/NubYStJ51/01rdyuzda0btUyoWhERKKL8kjqPOAGYBmw\nmmA2WQd2iTe00jNr3jd8uyh9qdVbLj8+oWhEROonyiOpK4A93f0/cQdTytydAYOfSqvToygRKSZR\n3pL6DPgh7kBK3S1DXksrb92yuR5FiUhRiXKHcTXwrpmNBZavqXT3S2KLqsQsXrKMcR9/nlb3f9f1\nSyYYEZH1FCVh3Ae8CXxM0Ich9XTmbx9JK193YS9NLCgiRSdKwmjk7gNij6REvTPh0xp1GqAnIsUo\nSh/GK2Z2npltZ2ZbrPkTe2QlwN2547HRaXWP33p2QtGIiGyYKHcYfcO/r06p02u1EVx+6zNp5V7d\nOtKkcaOEohER2TBRBu5p3qj1sHjJMj6fm/4mspZbFZFipvUwYlK9o/u2K05IJhARkRzRehgx+OeM\nuTXqdm69VQKRiIjkjrl73a1SdwjXw3D3nvGElPG8Xt9Yk3L8r+9NKw+9+Syabdw4oWhEpJyZGe6e\nk/f4tR5Gjj3/5kdp5b33aKNkISIlQeth5JC789jz76XV/e6CIxOKRkQkt7QeRg7d/kj6mIszjjlQ\nI7pFpGRkTBhm1h5otWY9jJT6g8ysibt/Fnt0RaSqahXvTUq/JL2775VQNCIiuZetD+NO4L+11P83\n3CYpBgxOf0p33YW9EopERCQe2RJGK3f/uHplWLdTbBEVoaXLVjDv60VpdZovSkRKTbaEsXmWbRvn\nOpBi1v9/Hk0r3zHwpIQiERGJT7aE8YGZnVu90szOASbEF1Jx+f6H5VStXJVWt+N2mptRREpPxoF7\nZtYKeA5YwboEsR/QGDjW3b/MS4Tr4inIgXvVB+ndc20/Wm3ZIqFoRETS5XLgXsa3pNz9K6CrmXUH\n9gyrX3L3N3Nx4lLw/Q/La9QpWYhIqar31CBJKcQ7jL6XP8iKqpVry/dedwpbb9E8wYhERNIlPTWI\nELwZlZosACULESlpShjr6Yrb0hdH+ss1/RKKREQkP5Qw1kNV1SrmVxt3se1W6rsQkdKmhLEebn1o\nVFr5j1dp3IWIlL7YE4aZ9TSzqWY23cyuqmV7PzP7KPzzjpn9JO6YNoS7M+GTWWl1bbfXuAsRKX2x\nJgwzawDcDfQAOgB9zWyPas1mAj91972AG4EH4oxpQw0d+X5a+YaLeycUiYhIfsV9h9EFmOHus9y9\nChgBHJ3awN3fd/c1HQLvAzvEHNMGqb5AUof22ycUiYhIfsWdMHYAZqeU55A9IZwDvBJrRBtg7OR/\np5XPOeHghCIREcm/KAso5UU4orw/kPGn8KBBg9Z+rqiooKKiIva4Ut0y5LW08hGH7JmhpYhIMior\nK6msrIzl2LGO9DazA4BB7t4zLA8E3N0HV2vXEXgG6JlpYaakR3rP/3oRF904fG35wL3bcXn/wxKL\nR0QkimIa6T0eaG9mbc2sMdAHGJnawMx2JEgWpxXyKn6X3fJ0WvnS0w5NKBIRkWTE+kjK3VeZ2UXA\nKILkNMTdp5jZ+cFmvx+4BtgCuMeCBbCr3L1LnHHV18qVq1i+omptedNmTWjYcKMEIxIRyT9NPhjB\n3Y+PYczYaWvLmsJcRIpFMT2SKgmpyQI0hbmIlCcljDp88tn8tPIlp6rvQkTKkxJGHa750/Np5W6d\nd0soEhGRZClhZLFseVVa+cfttksoEhGR5ClhZHHfk2+nlS/vf3hCkYiIJE8JI4u3P5iRVt6s+cYJ\nRSIikjwljAxmzv46rfxrDdQTkTKnhJHB/977clr5kH13TSgSEZHCoIRRC3fnv98vXVtu3aolwSB0\nEZHypYRRi8px09PKV57TI6FIREQKhxJGLe5+fExaeYdtNk8oEhGRwqGEUc2KqpVp5UP3r76irIhI\neVLCqGb4S+PTymcee2BCkYiIFBYljGpGjklfs3uTjZskFImISGFRwkiR+mYUwCm99k8oEhGRwqOE\nkWLoyLFp5d7dOyYUiYhI4VHCSPHm2KlpZa2qJyKyjhJGaNHi9MdRF/btllAkIiKFSQkjNOzF9MdR\n3bvsnlAkIiKFSQkj9Mb76Y+jGjTQpRERSaWfisCSpcvTyuef9NOEIhERKVxKGMDzb6SPvTh0fz2O\nEhGpTgkDeOb1iWllvR0lIlJT2SeMlStXpZX7HtUloUhERApb2SeMMeOmpZWPOKRDQpGIiBS2sk8Y\nDz37blpZc0eJiNSu7BNG6nTmWoZVRCSzsk4Y0z//Kq18Ys99E4pERKTwlXXCePLVD9LKWllPRCSz\nsk4YH06ZvfZzqy1bJBiJiEjhK9uEsXTZirTyab0PSCgSEZHiULYJ4/X3pqSVO+/ZNqFIRESKQ9km\njBEvp/dfaHS3iEh2ZZswlq+oWvv5sK4/SjASEZHiUJYJY+6ChWnlo7ppKVYRkbrEnjDMrKeZTTWz\n6WZ2VYY2fzKzGWY2ycz2jjum0e+m91+02bZl3KcUESl6sSYMM2sA3A30ADoAfc1sj2ptjgDaufuu\nwPnAvXHGBDByzLrpzJs1bRz36XKusrIy6RAKhq7FOroW6+haxCPuO4wuwAx3n+XuVcAI4OhqbY4G\nHgNw97HAZmbWKq6A3D2t3OfIznGdKjb6n2EdXYt1dC3W0bWIR8OYj78DMDulPIcgiWRrMzes+6pa\nO/5w3ysbHNCcr75LK3fdp90GH1NEpBzEnTByasIns3J+zJYtmuX8mCIipciqP6LJ6cHNDgAGuXvP\nsDwQcHcfnNLmXmCMuz8RlqcC3dz9q2rHii9QEZES5u6Wi+PEfYcxHmhvZm2B+UAfoG+1NiOBXwFP\nhAlmYfVkAbn7wiIisn5iTRjuvsrMLgJGEXSwD3H3KWZ2frDZ73f3l83sSDP7FFgC9I8zJhERWT+x\nPpISEZEcoYNBAAAG7ElEQVTSURQjvaMM/isVZtbazN40s3+Z2cdmdklY39LMRpnZNDN7zcw2S9nn\n6nDg4xQzOzy56ONhZg3MbKKZjQzLZXktzGwzM3sq/G7/MrP9y/ha/MbM/mlmk81smJk1LpdrYWZD\nzOwrM5ucUlfv725mncLrN93M7ox0cncv6D8ESe1ToC3QCJgE7JF0XDF+322BvcPPmwLTgD2AwcCV\nYf1VwM3h5x8DHxI8XtwpvFaW9PfI8TX5DfBXYGRYLstrATwC9A8/NwQ2K8drAWwPzAQah+UngDPK\n5VoABwN7A5NT6ur93YGxQOfw88tAj7rOXQx3GFEG/5UMd//S3SeFn78HpgCtCb7zo2GzR4Fjws+9\ngRHuvtLdPwdmUHOsS9Eys9bAkcCDKdVldy3MrAVwiLs/DBB+x0WU4bUIbQRsYmYNgY0Jxm+VxbVw\n93eA76pV1+u7m9m2QHN3Hx+2eyxln4yKIWHUNvhvh4RiySsz24ngN4n3gVYevj3m7l8C24TNMg18\nLBV3AFcAqZ1t5Xgtdgb+Y2YPh4/n7jezZpThtXD3ecDtwBcE32uRu4+mDK9Fim3q+d13IPhZukak\nn6vFkDDKkpltCjwN/Dq806j+dkLJv61gZkcBX4V3XNleqy75a0HwSKET8Bd370TwRuFAyvPfxeYE\nv1G3JXg8tYmZnUIZXossYvnuxZAw5gI7ppRbh3UlK7zNfhoY6u7Ph9VfrZljK7ydXBDWzwXapOxe\nStfnIKC3mc0EhgOHmtlQ4MsyvBZzgNnuvmblr2cIEkg5/rv4OTDT3b9191XAc0BXyvNarFHf775e\n16QYEsbawX9m1phg8N/IhGOK20PAJ+5+V0rdSODM8PMZwPMp9X3Ct0R2BtoD4/IVaJzc/bfuvqO7\n70Lw3/1Ndz8NeIHyuxZfAbPNbLew6mfAvyjDfxcEj6IOMLOmZmYE1+ITyutaGOl33fX67uFjq0Vm\n1iW8hqen7JNZ0j3+Ed8K6EnwttAMYGDS8cT8XQ8CVhG8DfYhMDH8/lsAo8PrMArYPGWfqwnefpgC\nHJ70d4jpunRj3VtSZXktgL0IfoGaBDxL8JZUuV6L68LvNZmgk7dRuVwL4HFgHrCcIHn2B1rW97sD\n+wIfhz9X74pybg3cExGRSIrhkZSIiBQAJQwREYlECUNERCJRwhARkUiUMEREJBIlDBERiUQJQwqa\nma02s1tTypeZ2bU5OvbDZnZcLo5Vx3lOMLNPzOyNavVtzaz6CpRRj/lObqITiU4JQwrdcuA4M9si\n6UBSmdlG9Wh+NnCOu/+sWv3OQL/1Ob67H1yP84vkhBKGFLqVwP3AgOobqt8hmNni8O9uZlZpZn8z\ns0/N7CYz62dmY83so3CKhDUOM7PxFizQdVS4fwMzuyVsP8nMzk057ttm9jzBtBzV4+kbLkgz2cxu\nCuuuIVi/YIiZDa62y03AweHss782szPM7PnwTmS0mW1iZqPN7IMw7t4ZvusYW7ew0tCUNjdbsMjQ\nJDO7pV5XXaQWsa7pLZIDDvwF+LiWH7i1tV2jI8HCUwsJFtt5wN33t2AFw4tZl4DauntnM2sPjDGz\ndgRz8SwM2zcG/mFmo8L2+wAd3P2L1BOb2XbAzeH2hcDrZtbb3f/XzA4FBrj7h9XiHQhc5u69w2Oc\nEe7/E3dfZGYNgGPc/Xsz25Jgmvs186ilfte9CRbK+TKMtSswNdx3j/DYLeq4diJ10h2GFDwPpnd/\nFPh1PXYb7+4L3H0F8BnB/DoQzJ2zU0q7J8NzfBq22wM4HDjdzD4kWJVsC2DXsP246ski1BkY48EM\nqquBYcBPU7Znm5491eseLIwEwf+fN5nZRwTzBG1vZtvUss84d5/vwTw/k8LvtwhYamYPmtmxwNKI\n5xfJSAlDisVdBH0Bm6TUrST8NxzOuNk4ZdvylM+rU8qrSb+zTv1N3cKyARe7+z7hn3YeLNADwToU\nmURNCtmkHv8UYCtgH3ffh2DK6qa17JP6XVcBDT2Y9rsLwTT5vYBXcxCblDklDCl0BuDu3xHcDZyd\nsu1zYL/w89EEM5bW14kWaEfQCT0NeA24MFyXBDPb1YLV7bIZB/zUzLYIO6z7ApV17LMYaJ5l+2bA\nAndfbWbdCRYMWiNrcgrj3dzdXyV4/NaxjlhE6qQ+DCl0qXcAtwO/Sql7AHg+fHT0Gpl/+882JfMX\nBD/smwPnu/sKM3uQ4LHOxPDOZQF1rHfs7l+a2UDWJYkX3f3FOs4/GVgdxv8INddpHga8ED6S+oBg\neuq6vtOa+hYE12bNHclvssUvEoWmNxcRkUj0SEpERCJRwhARkUiUMEREJBIlDBERiUQJQ0REIlHC\nEBGRSJQwREQkEiUMERGJ5P8Bzu3gnS9rJRgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7faa0a78a650>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cdf = Cdf(suite)\n",
    "thinkplot.Cdf(cdf)\n",
    "thinkplot.Config(xlabel='Number of trains',\n",
    "                 ylabel='Cumulative Probability',\n",
    "                 legend=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Cdf` also provides `Percentile`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(69, 869)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cdf.Percentile(5), cdf.Percentile(95)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** To write a likelihood function for the locomotive problem, we had\n",
    "to answer this question:  \"If the railroad has `N` locomotives, what\n",
    "is the probability that we see number 60?\"\n",
    "\n",
    "The answer depends on what sampling process we use when we observe the\n",
    "locomotive.  In this chapter, I resolved the ambiguity by specifying\n",
    "that there is only one train-operating company (or only one that we\n",
    "care about).\n",
    "\n",
    "But suppose instead that there are many companies with different\n",
    "numbers of trains.  And suppose that you are equally likely to see any\n",
    "train operated by any company.\n",
    "In that case, the likelihood function is different because you\n",
    "are more likely to see a train operated by a large company.\n",
    "\n",
    "As an exercise, implement the likelihood function for this variation\n",
    "of the locomotive problem, and compare the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# Suppose Company A has N trains and all other companies have M.\n",
    "# The chance that we would observe one of Company A's trains is $N/(N+M)$.\n",
    "# Given that we observe one of Company A's trains, the chance that we\n",
    "# observe number 60 is $1/N$ for $N \\ge 60$.\n",
    "\n",
    "# The product of these probabilities is $1/(N+M)$, which is just the\n",
    "# probability of observing any given train.\n",
    "\n",
    "# If N<<M, this converges to a constant, which means that all value of $N$\n",
    "# have the same likelihood, so we learn nothing about how many trains\n",
    "# Company A has.\n",
    "\n",
    "# If N>>M, this converges to $1/N$, which is what we saw in the previous\n",
    "# solution.\n",
    "\n",
    "# More generally, if M is unknown, we would need a prior distribution for\n",
    "# M, then we can do a two-dimensional update, and then extract the posterior\n",
    "# distribution for N.\n",
    "\n",
    "# We'll see how to do that soon."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Suppose I capture and tag 10 rock hyraxes.  Some time later, I capture another 10 hyraxes and find that two of them are already tagged.  How many hyraxes are there in this environment?\n",
    "\n",
    "As always with problems like this, we have to make some modeling assumptions. \n",
    "\n",
    "1) For simplicity, you can assume that the environment is reasonably isolated, so the number of hyraxes does not change between observations.\n",
    "\n",
    "2) And you can assume that each hyrax is equally likely to be captured during each phase of the experiment, regardless of whether it has been tagged.  In reality, it is possible that tagged animals would avoid traps in the future, or possible that the same behavior that got them caught the first time makes them more likely to be caught again.  But let's start simple.\n",
    "\n",
    "I suggest the following notation:\n",
    "\n",
    "*  `N`: total population of hyraxes \n",
    "*  `K`: number of hyraxes tagged in the first round \n",
    "*  `n`: number of hyraxes caught in the second round\n",
    "*  `k`: number of hyraxes in the second round that had been tagged\n",
    "\n",
    "So `N` is the hypothesis and `(K, n, k)` make up the data.  The probability of the data, given the hypothesis, is the probability of finding `k` tagged hyraxes out of `n` if (in the population) `K` out of `N` are tagged. \n",
    "\n",
    "If you are familiar with the hypergeometric distribution, you can use the hypergeometric PMF to compute the likelihood function.  Otherwise, you can figure it out using combinatorics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "from scipy.special import binom\n",
    "\n",
    "class Hyrax(Suite):\n",
    "    \"\"\"Represents hypotheses about how many hyraxes there are.\"\"\"\n",
    "\n",
    "    def Likelihood(self, data, hypo):\n",
    "        \"\"\"Computes the likelihood of the data under the hypothesis.\n",
    "\n",
    "        hypo: total population (N)\n",
    "        data: # tagged (K), # caught (n), # of caught who were tagged (k)\n",
    "        \"\"\"\n",
    "        N = hypo\n",
    "        K, n, k = data\n",
    "\n",
    "        if hypo < K + (n - k):\n",
    "            return 0\n",
    "\n",
    "        like = binom(N-K, n-k) / binom(N, n)\n",
    "        return like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0010248876709531896"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "hypos = range(1, 1000)\n",
    "suite = Hyrax(hypos)\n",
    "\n",
    "data = 10, 10, 2\n",
    "suite.Update(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEPCAYAAACQmrmQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcXGWd7/HPr/dOJ+nsW2cjK4R9NSpKCyKLShQdBWcG\nZbwjd5TR63jnouOMJHPvdWFmHEWuA8wgA44Kjo4aFQEXGhUhBEJIIPu+J52ks3V6rfrdP86p6tOV\nXqqX09Xd9X2/XvXqc556nnOeOoT61bOc55i7IyIiEoeCXFdARESGLwUZERGJjYKMiIjERkFGRERi\noyAjIiKxUZAREZHYxB5kzOx6M9tgZpvM7K5O8txrZpvNbLWZXdRdWTO70MyeN7NXzOxFM7ss7s8h\nIiI9F2uQMbMC4D7gOuBc4FYzOzsjzw3AXHefD9wB3J9F2XuAu939YuBu4B/i/BwiItI7cbdkrgA2\nu/tOd28BHgOWZORZAjwK4O4rgEozm9xN2SRQGW6PAfbG+zFERKQ3imI+fhWwO7K/hyB4dJenqpuy\nnwaeMrN/Agx4Uz/WWURE+slgHPi3LPL8BfApd59JEHC+FW+VRESkN+JuyewFZkb2p3Nm19ZeYEYH\neUq6KPthd/8UgLv/wMwe6ujkZqaF2UREesHds/nB3624WzIrgXlmNsvMSoBbgOUZeZYDtwGY2WLg\nmLsf7KTsT8Iye83sqrDMNcCmzirg7nq5c/fdd+e8DoPlpWuha6Fr0fWrP8XaknH3hJndCTxNENAe\ncvf1ZnZH8LY/6O5PmNmNZrYFqAdu76LshvDQfw7ca2aFQCPwsTg/h4iI9E7c3WW4+5PAwoy0BzL2\n78y2bJj+B0D3xoiIDHKDceBfYlBdXZ3rKgwauhZtdC3a6FrEw/q7/20wMTMfzp9PRCQOZoYPkYF/\nERHJYwoyIiISGwUZERGJjYKMiIjERkFGRERioyAjIiKxUZAREZHYKMiIiEhsFGRERCQ2CjIiIhIb\nBRkREYlN7KswD0XuTs2Lm3jptR2cv2A61125CLN+WcZHRCSvKMh04DcrNvDN7z0LwAtrtnOqoYn3\nv+OSHNdKRGToUXdZBnfnB0+tapf2X798hfqGphzVSERk6FKQybBj7xEOHT3ZLq2puYXfvLAxRzUS\nERm6FGQybNh+oMP051/dNsA1EREZ+hRkMqzf1hZk3v+OS0gN92/afoBjJ0/nplIiIkNU7EHGzK43\nsw1mtsnM7uokz71mttnMVpvZRd2VNbPHzGxV+NpuZqs6Om5vbN9dm96+/LzZLJwzBQAHXtu8r79O\nIyKSF2INMmZWANwHXAecC9xqZmdn5LkBmOvu84E7gPu7K+vut7j7Je5+CfBD4L/6o76JRJIDR9rG\nY6ZPGct586vS++u27O+P04iI5I24WzJXAJvdfae7twCPAUsy8iwBHgVw9xVApZlNzrIswAeA7/VH\nZQ8dPUkymQRgXGUFZaXFnDt3avr9dVvVkhER6Ym4g0wVsDuyvydMyyZPt2XN7C3AAXff2h+V3V97\nPL09dWIlAAtmT6agILhMuw/UcfxkQ3+cSkQkLwzGgf+e3Fp/K/3UigHYd+hYejsVZMpKi5k7Y0I6\nffOuQ/11OhGRYS/uO/73AjMj+9PDtMw8MzrIU9JVWTMrBG4GurwVf+nSpent6upqqqurO8175Fh9\nenvS+FHp7fmzJrF5ZxBctu6q5bJzZ3V1ShGRIaWmpoaamppYjh13kFkJzDOzWcB+4BaC1kfUcuAT\nwONmthg45u4HzexwN2WvBda7e5cDJdEg050jx9uCzPjKivT2vJmT0tvbIrPPRESGg8wf4MuWLeu3\nY8caZNw9YWZ3Ak8TdM095O7rzeyO4G1/0N2fMLMbzWwLUA/c3lXZyOE/SD92lQHURYLM2EiQmTNj\nYnp7q4KMiEjWYl8g092fBBZmpD2QsX9ntmUj793eX3VMqTvRdrPluEiQqZpUSWlJMU3NLdSdOM3R\n4/Xt3hcRkY4NxoH/nDl6PBpkRqS3CwoKmDO9bfB/yy61ZkREsqEgE2pobKapuQWA4qJCRpSVtHt/\nrrrMRER6TEEmdDTSVTZ29IgzHlI2JzKNede+IwNWLxGRoUxBJnTyVGN6u3JU+Rnvz5o2Lr29c9/R\nAamTiMhQpyATOlHfFmRGVZSd8X7VpLEUhK2bg0dO0NDYPGB1ExEZqhRkQqciQWbkiNIz3i8uLqRq\n8pj0/u4DdQNSLxGRoUxBJnTydNvjlUdXnNldBjBz2vj09o69GpcREemOgkzo5Km2hS9HVpzZkoH2\n4zK79mtcRkSkOwoyoWhLZtSIM8dkAGZFWjI7NcNMRKRbCjKh6JjMqJGdBJmpbS2ZHXuP4O6x10tE\nZChTkAm1m13WwcA/wISxI9M3aZ5ubG63arOIiJxJQSZ0Ktpd1sEUZgAzU5eZiEgPKMiETkfue6no\npCUDMHOqbsoUEcmWgkzodENbkCkvLe403+yqSEtmv1oyIiJdUZAB3J3Gppb0fldBJtqS2aWWjIhI\nlxRkgOaWVpLhTLHiokKKigo7zTtj6tj09t5Dx2htTcRePxGRoUpBBjjdGGnFZCzxn6mivJQJY0cC\nkEgk2Vd7PNa6iYgMZQoytB/0H1HWeVdZSrsuM935LyLSKQUZoCE66N9NSwba35SpcRkRkc7FHmTM\n7Hoz22Bmm8zsrk7y3Gtmm81stZldlE1ZM/tLM1tvZmvN7Mt9qWNDloP+KTO1hpmISFaK4jy4mRUA\n9wHXAPuAlWb2E3ffEMlzAzDX3eeb2RuA+4HFXZU1s2rg3cD57t5qZhPog/bdZVm0ZHRDpohIVuJu\nyVwBbHb3ne7eAjwGLMnIswR4FMDdVwCVZja5m7J/AXzZ3VvDcof7Usl205ezCDLTJo5JP8Ds0NGT\n7cqLiEibuINMFbA7sr8nTMsmT1dlFwBvNbMXzOwZM7usL5XsaUumuLiQaZPaHmCmLjMRkY7F2l3W\nS5ZFniJgrLsvNrPLge8DczrKuHTp0vR2dXU11dXVZ+SJBpnyLGaXQTAus+dg8HTMXfuPsmD25KzK\niYgMNjU1NdTU1MRy7LiDzF5gZmR/epiWmWdGB3lKuii7B/gvAHdfaWZJMxvv7mcMkESDTGd6OrsM\ngmnMf3hlK6CWjIgMbZk/wJctW9Zvx467u2wlMM/MZplZCXALsDwjz3LgNgAzWwwcc/eD3ZT9MXB1\nWGYBUNxRgMlWu5sxs5hdBhr8FxHJRqwtGXdPmNmdwNMEAe0hd19vZncEb/uD7v6Emd1oZluAeuD2\nrsqGh/4W8C0zWws0EQap3mpo6tmYDGTekFnXl9OLiAxbsY/JuPuTwMKMtAcy9u/MtmyY3gL8aX/V\nsSHSkinLsiUzefwoSkuKaWpu4cSpBo6dPM2YUSP6q0oiIsOC7vgnoyVTnl1LxsyYMaVtsUzd+S8i\nciYFGdqPyWSzdlmKHmAmItI1BRmgMTKFuaw0u5YMwCwtLyMi0iUFGaCppTW9XVqS/TCVVmMWEema\nggzQ1BwJMsXZB5noNObdB+rw8MFnIiISUJChfZApK80+yFSOKmf0yPLwGC0cOHyi3+smIjKU5X2Q\ncXdaIo9QLulBSwY0LiMi0pW8DzLRVkxxUSFm2Syd1kbjMiIinVOQae7doH9KtCWjacwiIu0pyPRy\nZllKtCWzWy0ZEZF2FGR6ObMsZcaUtiCz79AxWloSXeQWEckveR9kmiNBpqQk+7v9U8pKi5k8fjQA\nSXf2HtJimSIiKXkfZBqb25aU6U13GWh5GRGRzuR9kOlrdxloGrOISGcUZPo4uwxghqYxi4h0KO+D\nTHNLdEymty0ZPSVTRKQjeR9k2i0p08sgM3XCaAoLg0t55Fg99Q1N/VI3EZGhTkGmH7rLiooKqZo0\nJr2vwX8RkYCCTEvfB/4Bzpo+Ib29fc/hPtVJRGS4iD3ImNn1ZrbBzDaZ2V2d5LnXzDab2Wozu6i7\nsmZ2t5ntMbNV4ev63tav/X0yvQ8ycyJBZpuCjIgIAL3/Vs2CmRUA9wHXAPuAlWb2E3ffEMlzAzDX\n3eeb2RuA+4HFWZT9qrt/ta91bN9d1vObMVPmzJiY3t62u7ZPdRIRGS7ibslcAWx2953u3gI8BizJ\nyLMEeBTA3VcAlWY2OYuyPVsuuRNNLZGbMfvSXVbVNsNsz4G6drPWRETyVdxBpgrYHdnfE6Zlk6e7\nsneG3Wv/ZmaVva1gfwz8A5SXlTBtYlCNpDu7NPgvIhJvd1kvZdNC+Sbw9+7uZvZ/gK8CH+0o49Kl\nS9Pb1dXVVFdXt3u/v8ZkAM6aMZF9tceBYFxm3qxJfTqeiMhAqKmpoaamJpZjxx1k9gIzI/vTw7TM\nPDM6yFPSWVl3jw56/Cvw084qEA0yHenrUv9Rc6ZP4LlVWwDYtkfjMiIyNGT+AF+2bFm/HTvu7rKV\nwDwzm2VmJcAtwPKMPMuB2wDMbDFwzN0PdlXWzKZEyt8MvNbbCvbH2mUp7WaY7dYMMxGRWFsy7p4w\nszuBpwkC2kPuvt7M7gje9gfd/Qkzu9HMtgD1wO1dlQ0PfU841TkJ7ADu6G0d293xX9rH7rJIkNmx\n7witrQmKigr7dEwRkaEs9jEZd38SWJiR9kDG/p3Zlg3Tb+uv+kWDTEkfWzKjKsqYOHYUtXUnSSSS\n7DlYx+yqCd0XFBEZpvL+jv+W1rYnWRb3Q6tj7gx1mYmIpCjIRIJMX1syEMwwS9Gd/yKS7/I+yERv\nmuyPloyWlxERaZP3QaalNZne7pcg0667rJbWSEtJRCTf5H2QaU3075jMmFEjmDB2JBB0xe0+UNfn\nY4qIDFV5HWTcnUSirSVTVNQ/l2P+rMnp7U07DvbLMUVEhqK8DjLRQf+iokLM+mXNTeZHlpPZtPNQ\nvxxTRGQoUpAJ9UdXWcqCSJDZrJaMiOSxvA4yzS3xBJk5MyZQUBBc2r2HjlHf0NRvxxYRGUryOsi0\nb8n036UoLSlm1rRx6f0tu7RYpojkJwWZUH+2ZAAWRAb/N24/0K/HFhEZKvI6yLTGGGSig/9bdqol\nIyL5Ka+DTEt0TKYflpSJmj87OsPsIO7er8cXERkK8jvIxNiSqZo0hhFlJQCcrG/kwOET/Xp8EZGh\nIK+DTHNMA/8AZsaC2bopU0TyW14HmThbMgBnz2l7gOe6rfv7/fgiIoNdl0HGzP49sv3h2GszwFpi\nuk8mZdHcqentdVv29fvxRUQGu+5aMhdGtj8VZ0VyoTVjWZn+Nn/WJAoLg0u8r/Y4x0829Ps5REQG\ns+6CzLCeEtXfDyzLVFJcxLyZbbPM1GUmIvmmuyAz3czuNbNvRLbTr2xOYGbXm9kGM9tkZnd1kude\nM9tsZqvN7KJsy5rZZ8wsaWbjMt/LRlx3/EctiozLbNimmzJFJL909/P9ryPbL/X04GZWANwHXAPs\nA1aa2U/cfUMkzw3AXHefb2ZvAO4HFndX1symA9cCO3tar5S41i6LOmfuVH7069UArNumloyI5Jcu\ng4y7P9LH418BbHb3nQBm9hiwBNgQybMEeDQ83wozqzSzycBZ3ZT9Z4IguLy3lYt7dhkEM8yMoN9x\n++5aTjc0M6K8JJZziYgMNl0GGTPr8gvc3W/q5vhVwO7I/h6CwNNdnqquyprZTcBud1/bl2fADESQ\nqSgvZea08ezcdwQHNu44yMXnzIjlXCIig0133WVvJPii/x6wAuifp3p1rctzmFk58DcEXWXdllm6\ndGl6u7q6murq6vR+3LPLUhbNncrOfUcAWL91v4KMiAwqNTU11NTUxHLs7oLMFIIv81uBDwE/B77n\n7q9nefy9wMzI/vQwLTPPjA7ylHRSdi4wG3jVgmbMdOBlM7vC3c94DGU0yGSKe3ZZyqJ5U/nF714D\nYM2mPXzojMaciEjuZP4AX7ZsWb8du8spVe6ecPcn3f3DwGJgC1BjZndmefyVwDwzm2VmJcAtnDmG\nshy4DcDMFgPH3P1gZ2Xd/TV3n+Luc9z9LIJutIs7CjDdiQ78FxXGt/jB+fOr0k2tLTsP6SFmIpI3\nuv1mNbNSM7sZ+A/gE8C9wI+yObi7J4A7gaeB14HH3H29md1hZh8L8zwBbDezLcADwMe7KtvRaehl\nN177lkx83WWjKso4a8ZEIKjs61s0y0xE8kN3A/+PAucBTwDL3P21np7A3Z8EFmakPZCx32HLqKOy\nHeSZ09M6pbQm4h/4T7lwQRXbdgfPlVmzcQ9XnD871vOJiAwG3bVk/gSYT7CkzPNmdiJ8nTSzIb92\nfbvuspiDzPkLp6e312zcE+u5REQGi+7ukxnWqzS3DtDAP8A5c6ZQVFRIa2uCvYeOcbjuFBPGjoz1\nnCIiudbdKsxlZvY/zOw+M/uYmcX7TTzAomMycQ78QxDEzoksMbN2U+YkOxGR4ae7b9ZHgMuAtcCN\nwD/FXqMBNBA3Y0ZdsCDSZbZJXWYiMvx11zJZ5O7nA5jZQ8CL8Vdp4AzU7LKUCxZU8Z1w+9WNe3B3\n+rJigYjIYNddS6YlteHurTHXZcC1tCbT2wPRkpkzYwIjR5QCcPxkA9t2H479nCIiudTtQ8uiM8qA\nC4bT7LKWlra4GffsMoCCggIuWdS2iMFLr/d6AWkRkSGhuzv+C919dPga5e5Fke3RA1XJuAz0mAzA\nZefNTm+/rCAjIsPcsJ6i3J1cBJmLzp5OQUFw2bfurqXuxOkBOa+ISC4oyIQGKshUlJe2m8q8ap1a\nMyIyfOV5kGkb+B+I2WUpl547K7398uu7Buy8IiIDLc+DzMC3ZAAuO68tyKzesIeWyPI2IiLDSd4G\nmWQySTLZ1pIpjPmO/6hpEyuZMiGYN9HU3MLazbr7X0SGp7wNMpmtmIG8KdLMuDwyy+z51dsG7Nwi\nIgMpj4PMwN6ImemNF7U9oeDFtdvbLdYpIjJc5G2QaR7gGzEzLZg9mfFjKgA4dbqJ17fqQWYiMvzk\nbZBpTeRmZlmKmbH4wrbWzPOrtw54HURE4pa3QSbakslFdxnAGyNBZsWaHSQigU9EZDjI2yDTmqPp\ny1Fnz5nC2NEjADhxqoF16jITkWEm9iBjZteb2QYz22Rmd3WS514z22xmq83sou7Kmtnfm9mrZvaK\nmT1pZlM6Om5XcnWPTFRml9lzr2zJST1EROISa5AxswLgPuA64FzgVjM7OyPPDcBcd58P3AHcn0XZ\ne9z9Qne/GPg5cHdP69YcuQEyFwP/KW+6eG56+7lVW3VjpogMK3G3ZK4ANrv7TndvAR4DlmTkWQI8\nCuDuK4BKM5vcVVl3PxUpXwH0eDAj1wP/KefMmcKkcaMAON3YrOX/RWRYiTvIVAG7I/t7wrRs8nRZ\n1sz+j5ntAj4EfKGnFRsM3WUQdJm99fIF6f1nV27KWV1ERPpbd49fzoWsbr13978F/jYcq/lLYGlH\n+ZYubUuurq6muroaGByzy1Kuumw+P3jqZQBeXreL4ycbqBxVntM6iUj+qKmpoaamJpZjxx1k9gIz\nI/vTw7TMPDM6yFOSRVmA7wJPkEWQiYrOLsvlmAzAtEljWDB7Mpt2HCSZTPLcK1u48a3n57ROIpI/\noj/AAZYtW9Zvx467u2wlMM/MZplZCXALsDwjz3LgNgAzWwwcc/eDXZU1s3mR8u8B1ve0YoOluyyl\nOtJlVvOiusxEZHiINci4ewK4E3gaeB14zN3Xm9kdZvaxMM8TwHYz2wI8AHy8q7Lhob9sZmvMbDXw\nduBTPa1bS0t07bLc3y70povnpleC3rq7lu17Due4RiIifRf7mIy7PwkszEh7IGP/zmzLhunv72u9\noi2ZkuLcD02Nqihj8YVzeG5VcK/MU8+9zn//4FU5rpWISN/k/id8jgy27jKA6688N73925e2cLqh\nOYe1ERHpu7wNMs2tuV2FuSPnzJnCjCljgeBhZr99aXOOayQi0jd5G2RaWwZfS8bMeMebF6X3n3ru\nddw9hzUSEembvA0yuX5oWWeuunxBeoxo1/6jbNh2IMc1EhHpvTwOMtGWzOC5DBXlpbzl0rYZ2j+t\nWZPD2oiI9M3g+XYdYINtdlnUu6ovSG+/uGY7+2uP57A2IiK9pyDD4OouA5g5dRwXnxMsguDAT59R\na0ZEhqb8DTItg292WdRNb7swvf2bFRs4Wd+Yw9qIiPRO/gaZQTrwn3L+gipmV00AglbXL373Wo5r\nJCLSc3kcZAZvdxkE05mXXN02NvPzZ9fq5kwRGXIUZMjtQ8u68uaL5zFlwmgATp1u4gm1ZkRkiFGQ\nYXC2ZAAKCwt437WXpPd/+syras2IyJCSt0FmMD20rCtvvWw+k8erNSMiQ1PeBpnB9NCyrhQVFfL+\nd7S1Zpb/5lVOnW7KYY1ERLKXt0FmsM8ui3rrZfPTYzP1DU388OlVOa6RiEh28jjIDM5lZTpSVFTI\nre+8Ir3/89+u5dDRkzmskYhIdgb3t2uMBvOyMh1588VzmTdzEgCJRJLv/GxFjmskItI9BRkGf3cZ\nBPfNfPg9b0zv//7lLWzeeTCHNRIR6V7+BpnosjKFQ+MyLJo7lcvPm53ef+D7vyOZTHZeQEQkx2L/\ndjWz681sg5ltMrO7Oslzr5ltNrPVZnZRd2XN7B4zWx/m/6GZje5JnRKJJKlHgRnB/ShDxYff88Z0\ny2v7nsM89dy6HNdIRKRzsX67mlkBcB9wHXAucKuZnZ2R5wZgrrvPB+4A7s+i7NPAue5+EbAZ+FxP\n6tWSMX3ZzHr+4XJk6sRKbr724vT+d3/2InUnTuewRiIinYv7J/wVwGZ33+nuLcBjwJKMPEuARwHc\nfQVQaWaTuyrr7r9y91Q/0QvA9J5UaqiNx2R6zzUXMXViJQCnG5t55MfP57hGIiIdizvIVAG7I/t7\nwrRs8mRTFuDPgF/0pFJDbWZZppLiIv78j96S3v/dy5t5ce2O3FVIRKQTg/EbNuu+KzP7PNDi7t/t\nLM/SpUvT29XV1VRXVw/5lgzAhQunc+Wl8/j9y1sAuP/xZzlnzhRGVZTluGYiMtTU1NRQU1MTy7Hj\nDjJ7gZmR/elhWmaeGR3kKemqrJl9BLgRuLqrCkSDTEpzy9C5EbMr/+19V/Lapn0cO3ma4ycbePA/\nf8dnPnJtrqslIkNM6gd4yrJly/rt2HF/w64E5pnZLDMrAW4BlmfkWQ7cBmBmi4Fj7n6wq7Jmdj3w\n18BN7t7jhbyGyrpl3RlVUcZf3HpVev8Pr2xNt2xERAaDWIOMuyeAOwlmg70OPObu683sDjP7WJjn\nCWC7mW0BHgA+3lXZ8NDfAEYCvzSzVWb2zZ7Uazh0l6Vcdu4srn5D24S9f3n8WfbXHs9hjURE2sQ+\nJuPuTwILM9IeyNi/M9uyYfr8vtSpXZAZpA8s64nb3/smXt+yj4NHTtDY1MI/PvxLvvTp9wzJSQ0i\nMrwM3QGJPmg3u6xo6H8Rjygv4TMfuTZ9U+mOvYd5+Ed/yHGtREQUZIZ8d1nK3JkT+bP3vjm9//Rz\n63hmxcYc1khEJE+DzHCZXZbpuisX8aaL56b3/+XxZ1m/dX8OayQi+W74fMP2QLvZZcNgTCbFzPj4\nLVcxc+o4IFij7SsPPcXBIydyXDMRyVd5GWSGY3dZSnlZCZ/72A2MHlkOwMn6Rr74wC84Wd+Y45qJ\nSD7K+yAzHAb+M00aN4q7PnpdeiLAnoN1/N8HnqCxqSXHNRORfJP3QWa4tWRSzp4zhU/+8dXpNXo2\n7zzEPQ89RUtkPEpEJG55GWSG68B/pisvncdH339lev/VjXv42rd/3W5MSkQkTsP3G7YLLcN04L8j\nN7zlPD54w2Xp/Rde3cZXH/mVAo2IDIi8DDKJPOgui/qj6y7lXVddkN5fsWY79zz0NM2RR1CLiMQh\nL4NMU+TLdTgO/GcyMz7y3jdy09suTKe9vG4nX3zwF5xuaM5hzURkuMvLINP+oWXDvyUDQaC5bcli\n3nftJem0tZv28vmv/5jDdadyWDMRGc7yMshEB/7zaRFJM+ND77qCW995RTpt1/6jfO6ff8TOfUdy\nWDMRGa7yPsjkw5hMpve/4xLu/NDbKCgI/vMfPV7P5/75xzz3ytYc10xEhpu8DDKtw2yp/9542xsW\n8rf//UbKSosBaGpu4av//kse+fHzJBLJHNdORIaLvAwyTc2Rgf886i7LdOHC6Xzxf7yXKRNGp9OW\nP/Mqf/8vP+Po8foc1kxEhou8DDL5OPDfmVnTxnHP/3wfly6alU57bfM+Pv3l77NizfYc1kxEhoO8\nDDLR+0PycUwmU0V5KZ/72PV88IbL0svQnDrdxD0PPcU3v1ejNc9EpNfyMsioJXMmM+MD11/G3Z94\nN+PHVKTTf/3CBj71pcd56fWdOaydiAxVsQcZM7vezDaY2SYzu6uTPPea2WYzW21mF3VX1szeb2av\nmVnCzC7p6JhdabdAZh6PyXTk/AVVfPWuD7R7+NnhulN86cFf8I8P/5K6E6dzWDsRGWpiDTJmVgDc\nB1wHnAvcamZnZ+S5AZjr7vOBO4D7syi7Fngv8Gxv6tXc7o5/tWQyjRxRyl99+O186k+vZlRFWTr9\n+dVb+eT/fYzlz7yqtc9EJCtxt2SuADa7+053bwEeA5Zk5FkCPArg7iuASjOb3FVZd9/o7pshPYTQ\nI/l6M2ZPmBlvvWwB9/7NB7nq8gXp9NONzTzy4+f51Jce54VXt+HuOayliAx2cQeZKmB3ZH9PmJZN\nnmzK9kq+34zZE6NHlvPJP7maL3z8XUydWJlOP3D4BP/wraf5wjeWs27r/hzWUEQGs8H4M75XrZPO\nLF26NL1dXV1NdXW1Bv574cKF0/naZz/AL373Ot9/8iVONwYLa67bup+/u/cnnDd/GrfccDnnzJ2a\n45qKSE/V1NRQU1MTy7HjDjJ7gZmR/elhWmaeGR3kKcmibLeiQQYgkUiSTAZ3tBukH1Es3SsqKuTd\nb7uA6isW8P0nX+LJ369LX8vXNu/jbzcHweY911zMRWdPx6xffy+ISExSP8BTli1b1m/HjjvIrATm\nmdksYD9wC3BrRp7lwCeAx81sMXDM3Q+a2eEsykIPWz6ZM8v0RdhzoyrK+Oj7ruSdV13AD55+mWdf\n3EQyHJuej7aRAAASI0lEQVR5bfM+Xtu8jxlTxvLut13AWy6dr3EvkTxmcQ/cmtn1wNcJxn8ecvcv\nm9kdgLv7g2Ge+4DrgXrgdndf1VnZMP09wDeACcAxYLW739DBuT3z85041cDtn38ECGZRPfKl2/v/\nQ+eZ/bXH+eEvV7ULNimjR5bz9sVnc80bz2m3fI2IDF5mhrv3yy/w2INMLnUUZA7XneKOpf8BwNjR\nI/i3/31bLqo2LB08coInnn2NXz6/nqbmM1cJOH9BFW9/4zm84fyz8nZhUpGhoD+DTN71Y7Qf9M+7\njx+ryeNHc/vNb+IDN1zKr57fwM+fXcORY20Lba7dtJe1m/YyckQpb7p4LldeMo9Fc6eqy1JkGMu7\nb1mtWxa/ivJSllx9Ie+66nxeXreLX/1hPavW7STVpjx1uomnn1vH08+tY/yYCt588TyuvGQec2ZM\nUMARGWbyLshEl/kvLcm7jz+gCgsLuOL82Vxx/mwO153i1y9s4DcrNrR73PORY/Usf+ZVlj/zKuPH\nVHDZubO5/PzZnDdvmrrURIaBvBuTWbNxD8u++TMAzps/jWV33pSLquUtd2fd1v38ftUW/vDKVk6d\nbuowX2lJMRefM4NLF83kgoXTmTB25ADXVCR/aUymDxojLZmykuIc1iQ/mRnnzpvGufOm8dGb38ya\nTXv53cubefn1XdQ3tAWcpuYWXnh1Gy+8ug2AaRMruWDhdC5YOJ3z5k+jorw0Vx9BRHog/4JMU3N6\nu7RUQSaXiooKuWTRTC5ZNJPW1gQbth9g5dqdrHxtBwePnGiXd1/tcfbVHufJ37+OAWfNmMjZZ03m\n7DlTWTh7slo6IoNUHgaZtpZMuYLMoFFUVMh586s4b34VH3nvG9lz8Bgr1+5g7aa9rN+2v92sQAe2\n7a5l2+5anvjtawBMGDuShWdN4Zw5U5g3cyKzpo3X7EGRQSDv/i9siDzlUd1lg5OZMWPKWGZMGcvN\n115Mc0srG7cfZM3GPby6cQ/bdteSOZJ4uO4Uh+u28NyqLQAUmDF9yljmzJjI3BkTmDN9IrOrxlOm\nHxYiAyrvgkz0UcJlZfrCGQpKios4f0EV5y+o4o/f/QbqG5rYuP0gG7YdYOOOA2zacajd1HSApDu7\n9h9l1/6j1Ly4EQjWH5o6sZKZU8cxfeq4MJCNo2pSJUWazi4Si/wOMprCPCRVlJemx3IAWlsT7Nh7\nhA3bD7Bp5yG2765lX+3xM8o5bWM7rNmeTi8oKGDaxEpmTB1H1eQxTJtYydSJlUyZMJpRFWW6d0ek\nD/LuW7YxstxJeWlJDmsi/aWoqJB5syYxb9akdNrphma27z3Mtt2H2banlm27D7P3YN0Z3WwAyWSS\nPQfr2HOw7oz3RpSVBAEnDDzTJlYyefxoJowdybjKERQUaBVvka7kX5CJtmRK8+7j540R5SXpqdIp\njU0t7DlQx+4Ddew+cDT4u7+O2rqTnR7ndGMzW3fXsnV37RnvFRQUML6ygonjRjJh7EgmjRvFhLEj\nmRj+nTBmpMaAJO/l3bdsQ2NbkCnVwH9eKSstPqPFA9DQ2Myeg3Xs2n+U/YeOs//wCfbXHufA4RMd\nLvSZkkwmqa072WWQKi8rYeyocsaNqWDM6BGMG13B2Mq2v2NGj2B8ZYWCkQxbeRdk2nWXaeBfCALB\n/FmTmT9rcrt0d6fuxGkOHD7B/tpjHKg9wb7a49QePUlt3SlOnGro9tgNjc00NDZ3OEYUVVpSTOXI\nMkaPLGd06m9FGaNHllE5spxRqb8Vwd/ysmKNFcmQkH9BRlOYJUtmxrjKCsZVVrCog8dKN7e0Ult3\nitqjJzkc/j2U3j7F0RP1JBLJrM7V1NzCoaMtHDraeasoqrCwgNEVZYwcUUrFiFJGlod/R5RSUV5K\nRXlJh++NHFGq+4dkQOXdv7ZokFFLRvqipLiIqkljqJo0psP33Z2T9Y3UnTgdvI6f5uiJeuqOn6bu\neD11Jxs4eqyeupOnaY3cbJqNRCKZPm5PFRUVUlFeQnlpMeVlwd8RZSWUlQV/O0ovLy1hRFkxZeHf\n1Pt6fLl0J++CTH1D27Iy5WWaXSbxMbOw+6ucWdPGd5rP3alvaObEqQZOnGrkRH1j2/apBo6fauBk\nfSPHTzVy8lQjx081nHFfUE+0tiY4frKB4ye77+7rTmFhAaXFRZSVFlNaUkRpSfi3uCj4W1rU9n5x\nEaWpvyVFlJUWURJ9r6SI4uIiSooLKS4qpCSyra7BoSuvgoy7czzSjz66oiyHtREJmFm6K2vapO7z\nQ9C9duJUI/UNTZw6HbzqG5qob2imPtw/1dCU3q4/3cSphmbqG5qy7sLLRiKR5HSimdONzd1n7oPC\nwgJKiovC4FNISVHhGQGpuKiQ4uLUdpA/la+4qLDtVVxAUWEhhYUFFBUWUBSmF6X2CwspKgrSM9NS\n+TR1PXuxBxkzux74GlAAPOTuX+kgz73ADUA98BF3X91VWTMbCzwOzAJ2AB9w965HVgmmo6b+Byst\nKdaMHhmySkuKmTiumImM6lE5d6e5pZX6hmYamlpoSP1tagknKbRwurGZxqbgb0fpDU3NnG5sobGx\nucP7juKQSCRpSDTT97ZX/ygwo7AwDDrtglEYnCJpqWBWYAUUFhqFhYUUFFjwXkGYVhDkKywooLAg\nOHZhGMwKC4yiaJlIueD9VJpRVFSYrltwTqOoKDxOWK6gwNLnSOVNpaf2+1OsQcbMCoD7gGuAfcBK\nM/uJu2+I5LkBmOvu883sDcD9wOJuyn4W+JW732NmdwGfC9O6FO0eyLdWTE1NDdXV1bmuxqCQz9fC\nzMIureAHVl+uhbvT0pqgqbmVpuZWGptbaG5upTG139RCc0trB++30NjUSlNLK82R95qaW2lpaaWl\nNUlzaystLQmaWxM9Hq/qrdo9m5g4fUFWeZPuJFsT7RZulY7F3ZK5Atjs7jsBzOwxYAmwIZJnCfAo\ngLuvMLNKM5sMnNVF2SXAVWH5R4AasggytZEnMo4bU9GXzzXk5PMXayZdizZ9uRZmFo6bFDEqxv+d\nUsGspTVBc0vqbxiEWlrbpbW2JtuntQbbran3E0laW5O0JhK0tCZJJBK0JpK0tCb4zcYa5s+6sl16\nW95wP5EcsKA3XMQdZKqA3ZH9PQSBp7s8Vd2UnezuBwHc/YCZddqT/cUHfpHefnndzvT29MkdzwgS\nkcElGswqymM8Ue3LLP2rm7vN5u4kk05rJBC1BaEEiUSSlpa24JV0JxEGqEQiSSKZJJlMkkg4iWRb\nWmq/NRG+n3QSrQkSybbzRfOltpPJoA5JD/aDfGF6eM5kMjx2mC/YD/KkthPhdjLZf2N2MDgH/nsz\njaTTruFoYImaM31iL04jIvnOzMKxlQKG4/NZ3Z2Ce/+ifw8Y1wtYDDwZ2f8scFdGnvuBD0b2NwCT\nuyoLrCdozQBMAdZ3cn7XSy+99NKr56/+igNxt2RWAvPMbBawH7gFuDUjz3LgE8DjZrYYOObuB83s\ncBdllwMfAb4CfBj4SUcnd3dNrhcRyaFYg4y7J8zsTuBp2qYhrzezO4K3/UF3f8LMbjSzLQRTmG/v\nqmx46K8A3zezPwN2Ah+I83OIiEjvWNitJCIi0u+G5W2rZna9mW0ws03hfTTDmplNN7PfmNnrZrbW\nzD4Zpo81s6fNbKOZPWVmlZEynzOzzWa23szekbvax8PMCsxslZktD/fz8lqEtwT8Z/jZXjezN+Tx\ntfi0mb1mZmvM7DtmVpIv18LMHjKzg2a2JpLW489uZpeE12+TmX0tq5PHOfCfixdB4NxCsBpAMbAa\nODvX9Yr5M08BLgq3RwIbgbMJuhX/V5h+F/DlcHsR8ApBd+ns8HpZrj9HP1+TTwP/ASwP9/PyWgD/\nDtwebhcBlfl4LYBpwDagJNx/nGA8Ny+uBXAlcBGwJpLW488OrAAuD7efAK7r7tzDsSWTvgHU3VuA\n1E2cw5a7H/BwKR53P0Uw+246wed+JMz2CPCecPsm4DF3b3X3HcBmzrx/acgys+nAjcC/RZLz7lqY\n2WjgLe7+MED4GY+Th9ciVAhUmFkRUA7sJU+uhbv/Hsh8vniPPruZTQFGufvKMN+jkTKdGo5BprOb\nO/OCmc0m+MXyAhk3rQKpm1Yzr9Fehtc1+mfgrwmmYqbk47U4CzhsZg+HXYcPmtkI8vBauPs+4J+A\nXQSf67i7/4o8vBYRk3r42asIvk9TsvpuHY5BJm+Z2UjgB8CnwhZN5qyOYT/Lw8zeCRwMW3ZdTWEf\n9teCoLvjEuD/ufslBLM3P0t+/rsYQ/DLfRZB11mFmf0xeXgtuhDLZx+OQWYvMDOyPz1MG9bCLoAf\nAN9299R9QwfDdeAIm7qHwvS9wIxI8eF0jd4M3GRm24DvAVeb2beBA3l4LfYAu939pXD/hwRBJx//\nXbwd2ObuR909AfwIeBP5eS1SevrZe3VNhmOQSd8AamYlBDdxLs9xnQbCt4B17v71SFrqplVof9Pq\ncuCWcHbNWcA84MWBqmic3P1v3H2mu88h+G//G3f/U+Cn5N+1OAjsNrPU0sLXAK+Th/8uCLrJFptZ\nmZkZwbVYR35dC6N9675Hnz3sUjtuZleE1/A2OrkRvp1cz3qIaSbF9QQzrDYDn811fQbg874ZSBDM\npHsFWBVeg3HAr8Jr8TQwJlLmcwSzRtYD78j1Z4jpulxF2+yyvLwWwIUEP7xWA/9FMLssX6/F3eHn\nWkMw0F2cL9cC+C7BI1OaCALu7cDYnn524FJgbfjd+vVszq2bMUVEJDbDsbtMREQGCQUZERGJjYKM\niIjERkFGRERioyAjIiKxUZAREZHYKMjIkGZmSTP7h8j+Z8zsC/107IfN7Ob+OFY353m/ma0zs19n\npF9lZj+N+/wicVKQkaGuCbjZzMbluiJRZlbYg+wfBf6bu1/TwXu9vpHNzPT/t+Sc/hHKUNcKPAj8\nVeYbmS0RMzsZ/r3KzGrM7MdmtsXMvmRmHzKzFWb2ariURsq1ZrbSgofgvTMsX2Bm94T5V5vZn0eO\n+1sz+wnB8i2Z9bk1fODTGjP7Upj2dwTP+njIzL7SwecbZW0PHft2WOZtZvajyHHfbmY/TH1GM/tH\nM3uFYBmVvzOzF8Nz3h/mKQzT3hruf8nM/ne4fUl4bVaa2S8ia1t90oKHnq02s+9m959GhOG5rIxe\n+fMCThA8qG07MAr4DPCF8L2HgZujecO/VwFHCZY2LyFYSPLu8L1PAl+NlH8i3J5HsPx5CfDnwN+E\n6SUEy7bMCo97EpjZQT2nAjsJljEpAH4N3BS+9wxwcQdlriJ4BshUgjWn/gC8KXxvHTA+3P4OcGO4\nnQTeFzlGdKmQR4F3htuLCALhNcDLBM9aKQKeixz3A8BD4fZeoDjcHp3r/+56DZ2XWjIy5HnwWINH\ngE/1oNhKdz/k7s3AVoK1myBYl2l2JN/3w3NsCfOdDbwDuC1sLawgCBzzw/wvuvuuDs53OfCMB6sA\nJwkCw1sj73f2WIIX3X2/uzvB+mOpun0b+JPwkbmLgSfD9FaCNcpSrjGzFyx47O7bgHPDz7OO4Mmh\nPyN4cmYCWAicB/wy/GyfJ1gWH+BV4Lvh8viJTuoqcoaiXFdApJ98nWBh0Icjaa2EXcLhqrElkfea\nItvJyH6S9v9fRMdELNw34C/d/ZfRCpjZVQTPbOlMV8+36Uy0nolI3f6dYGXpJuA/w8AF0BgGJMys\nFPh/wCXuvs/M7gbKIsc7n6ClNDlSv9fc/c0d1OOdBEHxJuDzZnZe5JwinVJLRoY6A3D3OoJWx0cj\n7+0ALgu3lxCsuttTf2SBuQRPmtwIPAV8PHyGD2Y2P3ziZFdeBN5qZuPCSQG3AjW9qA8A7r6fYFXd\nz9M+sEYDWRlBUDxiwQPt3p/OFIxVjSUIHPdZ8KjmjcBEM1sc5ikys0VhkZnu/izBQ89GE3RRinRL\nLRkZ6qItjX8CPhFJ+1fgJ2HXz1N03sroagbXLoIAMQq4w92bzezfCLqtVoUtpEN086xzdz9gZp+l\nLbD8zN1/lsX5u6rnd4AJ7r6xozzuftzM/pVg7GV/+Dkws/HAF4GrwxbONwiWbb/dzP4IuDfshisE\nvmZmm4D/CAORhXlPZFlnyXNa6l9kiAqDwyp3f7jbzCI5oiAjMgSZ2UvAKeBad2/JdX1EOqMgIyIi\nsdHAv4iIxEZBRkREYqMgIyIisVGQERGR2CjIiIhIbBRkREQkNv8flpjUPlGmFpcAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7faa0a6f6290>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "thinkplot.Pdf(suite)\n",
    "thinkplot.Config(xlabel='Number of hyraxes', ylabel='PMF', legend=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior mean 185.570957948\n",
      "Maximum a posteriori estimate 50\n",
      "90% credible interval (36, 618)\n"
     ]
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "print('Posterior mean', suite.Mean())\n",
    "print('Maximum a posteriori estimate', suite.MaximumLikelihood())\n",
    "print('90% credible interval', suite.CredibleInterval(90))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "from scipy import stats\n",
    "\n",
    "class Hyrax2(Suite):\n",
    "    \"\"\"Represents hypotheses about how many hyraxes there are.\"\"\"\n",
    "\n",
    "    def Likelihood(self, data, hypo):\n",
    "        \"\"\"Computes the likelihood of the data under the hypothesis.\n",
    "\n",
    "        hypo: total population (N)\n",
    "        data: # tagged (K), # caught (n), # of caught who were tagged (k)\n",
    "        \"\"\"\n",
    "        N = hypo\n",
    "        K, n, k = data\n",
    "\n",
    "        if hypo < K + (n - k):\n",
    "            return 0\n",
    "\n",
    "        like = stats.hypergeom.pmf(k, N, K, n)\n",
    "        return like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.046119945192892968"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "hypos = range(1, 1000)\n",
    "suite = Hyrax2(hypos)\n",
    "\n",
    "data = 10, 10, 2\n",
    "suite.Update(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior mean 185.570957948\n",
      "Maximum a posteriori estimate 50\n",
      "90% credible interval (36, 618)\n"
     ]
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "print('Posterior mean', suite.Mean())\n",
    "print('Maximum a posteriori estimate', suite.MaximumLikelihood())\n",
    "print('90% credible interval', suite.CredibleInterval(90))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
